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PREFACE 

Background.  This  work  is  a  continuation  of  a  previous  work  effort  (1)  with 
the  objective  of  developing  a  viscoplastic  constitutive  model  for  soils  and 
rocks.  In  the  previous  work,  the  inviscid-plastlc  cap  model  of  Sandler  and 
Rubin  (CAP75)  was  reformulated  into  a  Perzyna-type  elastic/viscoplastic 
model  (2).  In  addition  to  the  theoretical  development  a  numerical  solution 
algorithm  was  developed  to  compute  six  dimensional  stress  histories  from  an 
arbitrary  strain  loading  schedule  and  vice  versa.  The  algorithm  was 
embodied  in  the  computer  program  "VPDRVR"  which  employs  a  Crank-Nicolson 
time  integration  scheme  and  a  Newt on -Raphson  iterative  solution  procedure. 
Numerical  studies  were  performed  to  validate  the  program  and  assess  the 
accuracy  for  various  options  of  the  time  integration  scheme.  The  effect  of 
the  model  fluidity  parameters  was  illustrated  for  triaxial  stress  and 
uniaxial  strain  loading  for  a  well -studied  sand  material  (McCormick  Ranch 
Sand).  Lastly,  a  finite  element  solution  methodology  incorporating  the 
viscoplastic  model  was  presented.  It  was  concluded  that  the  elastic- 
viscoplastic  model  shows  great  promise  for  capturing  the  viscoplastic 
nature  of  many  geological  materials.  Recommendations  for  future  advance¬ 
ment  of  the  model  were  to  incorporate  a  viscoplastic  tension-cutoff  cri¬ 
terion  and  to  establish  parameter  identification  techniques  with 
experimental  data.  Herein  lies  the  Impetus  of  this  study. 

Objective.  As  indicated  above,  this  report  addresses  two  main  areas:  (1) 
formulation  of  a  viscoplastic  tension-cutoff  model  to  be  Incorporated  into 
the  viscoplastic  cap  model,  and  (2)  development  of  parameter  Identification 
procedures  and  guidelines  for  the  cap  viscoplastic  model. 

Scope  and  Approach.  Part  I  of  this  report  deals  with  tension  cutoff.  The 
underlying  motivation  for  introducing  tension  cutoff  stems  from  the 
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recognition  that  soils  and  rocks  usually  exhibit  abrupt  changes  In  the>r 
stress-strain  behavior  when  loaded  in  tension,  i.e.,  rapid  tensile  stress 
release  as  micro-cracking  or  particle  separation  occurs.  To  this  end,  the 
Jj  (first  stress  invarlent)  tension-cutoff  criterion  proposed  by  Sandler 
and  Rubin  (3)  is  adopted  for  this  study  and  recast  Into  viscoplastic  for¬ 
mulation.  Here,  separate  fluidity  parameters  are  assigned  to  the  tension 
cutoff  domain  to  permit  independent  control  on  the  rate  of  tensile-stress 
release. 

For  the  sake  of  completeness.  Part  I  presents  a  concise  review  of  the 
viscoplastic  cap  model  prior  to  Introducing  the  viscoplastic  tension 
cutoff.  Following  the  theoretical  tension-cutoff  development,  a  numerical 
solution  algorithm  is  presented  for  the  entire  viscoplastic  model  Including 
tension  cutoff.  This  algorithm,  which  predicts  stress  histories  from 
strain  loading  and  vice  versa,  is  an  extension  of  the  previous  "VPORVR" 
program.  Input  instructions  for  the  new  "VPDRVR"  program  with  tension 
cutoff  is  given  in  the  Appendix.  Part  I  concludes  with  an  illustrative 
example  comparing  the  tension-cutoff  numerical  solution  with  an  exact  solu¬ 
tion.  Also  presented  is  a  critique  of  the  tension-cutoff  criterion  and 
recommendations  for  future  enhancements. 

Part  II  of  this  report  addresses  the  parameter  identification  problem. 
We  begin  by  Illustrating  the  Influence  of  various  model  parameters  on  the 
model's  performance.  With  these  Insights,  a  set  of  guidelines  are 
established  to  aid  In  parameter  Identification.  For  Identification  pur¬ 
poses,  experimental  data  is  classified  into  two  categories  "Ideal"  and 
"non-ideal".  The  former  implies  a  well -designed  set  of  experiments  espe-  ' 
daily  conceived  to  ease  the  identification  problem.  A  hypothetical 
example  of  an  Ideal  experiment  Is  presented  along  with  a  step-by-step  Iden- 
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tificatlon  procedure.  The  ’non- ideal  *  taper taeats.  «*tc*  IMl|  W»» 
applies  to  nost  eiisting  data.  Mi  «at  lend  itself  la  a  *d«#act*  ««••* 
tl  Heat  ion  procedure.  Hara  a  trial-and-errar  approach  asiap  »e  #tt* 
propraa  Is  suggested  and  threa  neap  I  at  of  actual  eaper inatal  data  ant 
investigated. 
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Although  more  elaborate  functional  forms  for  <|>  may  be  established  (2), 
the  forms  given  by  Equations  9  and  10  appear  to  suffice  for  many  geological 
materials  (4,5). 

Specification  of  the  plasticity  yield  function  f  Is  patterned  after  the 
Invlscld  cap  model  (3)  wherein  Jj,  the  first  stress  Invariant,  and  J'2,  the 
second  deviator  stress  Invariant,  are  used  to  define  the  current  static 
yield  surface,  as  Illustrated  In  Figure  1.  Here,  the  static  yield  surface 
Is  divided  into  three  regions  along  the  Jj  axis;  the  failure  surface  region 
(T  >  Ji  >  L) ,  the  cap  surface  region  (Jj  i  L),  and  the  tension-cutoff 
region  (Jj  >  T). 

Failure  Surface.  The  failure  surface  is  a  non-hardening,  modified 
Orucker-Prager  form  with  a  yield  function  defined  by: 

♦F(Jr  J‘2)  -  -  (A  -  C  exp(BJj))  (11) 

where  A,  B,  and  C  are  plastic  material  constants  (A  >  C).  This  yield  func¬ 
tion  is  used  to  define  viscoplastic  flow  (Equation  3)  whenever  Jj  Is  In  the 
range  T  »  jj  >  i.  The  failure  surface  forms  a  boundary  along  which  the  cap 
surface  can  move  (harden/soften). 

An  alternative  form  for  fF  Is  the  standard  Orucker-Prager  surface  given 
by: 

fr(Jr  J’2)  *  -  (A  -  BJj)  (12) 

where  A  and  •  are  material  constants.  The  first  form  is  generally  preferable. 
Cap  Surface.  The  cap  surface  Is  a  hardening  surface  in  the  shape  of  an 
ellipse  quadrant  when  plotted  in  space  (Figure  1).  It  is  defined 

In  a  "squared*  form  with  the  noraelltlng  constant  f0  (stress  unit)  as: 
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Figure  i.  i i lustration  of  steady-stele  pUstHH?  *«•'*»<?* 
For  cap  model . 


8 


»cUr  J  j.  *1  *  W  ^  •  VU  '  .  ...  .  t  ♦’  *(  Uli 

rn#  mmmIi  mi  m>wu»>  v  «n»  »  *>«  iw»hi«h  •»  u*  ^  *m*  #u* 

IMIIM  IN*  twM«M  Ml  W»A4*  i*4  -  ft  »*  *!.*•“,  N**-  *M  UaaNPM** 
t«l*«  )  »%  1*4*%  l»*  «A#  W *«K«  M4KVM4*  tM  •««>«**  wrf *t« . 

I  aa*  I  #'•  Mi*l*l  »y  Hf  «#'•*»•*  *A*lMrte4  t  «Mt*  ««Mw%  <*• 

C*er«  •  '  tn«  At  :»*-*•  'a**9  »«>a9N  Mw  ««*  *»**+i*. 

l  «  -  4  i«3 

*  •  •  C  «40(*t)  **  t*»  *Mt*«a*  «“»|W  »•!'**  #Ml  V  4  H  tM 

*  »*  «|4«<,><4  4*  4  inMUH  ,  %* 

Ip.  Mm***'.  %«•»*»«  t|4<',‘»4  ft.  (*  «  .  M*  N»»MU  *) 

C#f  *f  tn  *4  I«M'I<1||  •» 

»t»»int»ry  *  ♦  M«l  <•»«*  «  •*  •»■  M(«a|t|M«»  *» 

(MlfMttM.  ‘<  It  «%  f«««*  ty 

«{<}  *  !*(*/»«  ft/O  [}%) 

•Mft  t  iM  0  iMMit*  «M«M>  UN*  »aM»N  IN 

CAP  ul  *I,J.  tM  ft*'**!**  #4tWMt  «  »t  A*  MCMwUfHM  #f 

t«tM4t*U  |<«t»  if 

* 

* 

*  *  «o  *  j0  «  41  (?*) 

#NM  ‘  *•♦•(«»  *  *„  *  *  °V  (l») 

“  t  li  li 

tMf  tally  »M  («  H  )«c«tH  N  N*tlf|r*N  t  *  t*.  «U*  • 

Mil  M*9»Uv*  <*•!«•  «N»  IN  J)  •ftft*.  M»  l«  l««  r«*<m  tM  MHMl 
mImi  of  l  (ttaattaa  l«)  aM  (laaat IS).  UM>  IN  cap. 

l.t.,  J|  <  l  0*4  fc  '  0.  *<U«|Uttk  tin  KCKI  «*<ck  Mflttlfly 
incfMiMt  .  I.  AM  l;  tMraty  >*»— tN  cap  at  ft  amt  fa  tM 


* 


d<  fm<l  im. 

Atilt  1 1  f  )  44$  #  » '  <V  *i#  *mi  *  '  »  *4K  V«*  *4©  4  ‘v»  1  t*e  vlMl*  |W»M» 

EM  C4|>  It  fulfill  t***  !*«i>S  «*  '•*«  ’  »f*  Vw*»*i*  W 

in  jrdwr  Id  limit  «>  •».«•*#-  !%s>  »»U*U  »•»  H  »•»*  *M«  It 

itinwn#*.  u  naehmm*  ♦*♦**#  in**  »!♦•**  j*  *  l  «*»•>$ .  tews  *mr  sm'.»  (net 
re**>).  -.  iwt#  tee*  w  i»«.#v*wr*s  •  *•'•«**■*  WOw*  •  i  {C*i(«»n*  tw 

Iwn-jtuti  *, ■»?!**.«  ld©‘J>n»jj  J'»'i  Is  l*«  »« 1  *t  *  w»  ** *4 

-  *4 «  •  -  *  •■  .  i  U*) 

*V-  •  *0  >  .  *»  t  V 

*  ■  *  V  /  / 

*h*n  (am  in«f*me*t  it  KtuMkitln*  **  I  «»««l  »•*  *•<**»%  less  nape- 

t've.  r»|rt((l*<)  |Kf  (tt,  tMlHi,  s  >»  **  >MM<*t  ©♦«**  v*ch 

cn«e  in#  (w  tt'Mwur  v  «m  m  i*i«m  r*«i*'  i*©»  in  (t»w  *#)«♦  •» 

Jt,  HtA<*  l^t  wtiiM  *«i«f  •»  «  •***  >ot<in|  (fc»  U'U*t  «I  Vwi‘ 

fee#  f % 

*c  -  •  U  111) 

where  *t  l|  *(tn»w«  rtm  I  ©Ml  'on  |(  I If  IfUitf  l  »  Ji  l  t  .  (  *  J,  •  A/J*. 

*  '•  t  1  ? 

This  spec*#'  nerdwn'np  r*»t#  for  toll*  neips  to  limit  tutultf  dlle- 
1 9fKf  (I.*..  *©l*met'it  •♦p©n*»04»  '*'©'♦*11/  associated  wit*  1©4(3'«5  th* 

ft'lu't  |A(|  UM<«*  b)  a’,  iowtn©  {b»  cap  to  bf  Kll«|((<  loan*' 

when  SwOseouent  compression  toedi#©  *s  encownt er*d. 

Another  restriction  imposed  on  •$  ("tl  tt  •»  *♦»♦<■  iMatea  to  brcanr 
9reeter  then  its  initial  value  . 

The  for*9©ln9  completes  the  viscoplasttc  cep  model  review  and  sets  the 
st©9e  for  fonewletin9  the  viscoplestic  tension-cutoff  model. 

Tension  Cutoff  Theoretical  Development. 

In  this  section  we  present  e  viscoplastic  tension-cutoff  model 
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(W  MM  IIMH*  <rHe#lO*  K»«M  Kr  V*MlK  Ml  1^11 

for  too  UNiuif  coo  w*i  (}).  Altoemgp  tooir  (Mikioi  dee*  not  oeefera 
to  any  rigorous  froc t»re  IMory.  it  «t  to  Ot  Mf«ruu  for  (At 

purpose  of  rfftiljr  reduces  too  stiffens*  of  those  oloooot*  eipononcmg 
tension.  further  cringed  o«0  niwUM  of  too  tension  cutoff  Mil  in 
fl«oi  M  too  (onl«ft*t  portion  of  ims  development. 

Tension  Cotoff.  TOo  toot ion-cutoff  criterion  employed  *«  too  i*tiu<4  cop 
00O0I  tf  triggered  0O00  Jj  >  !.  Moot,  I  H  l  materiel  constant  npnw*( ■ 
top  too  threshold  of  nlooiru  teamen  tinti  ot  Uhich  Hoot  stress 
roliiiM  occur  Poo  to  tension  iooofo.  ipeof  Kelly,  OoMKf  o  ttnit 
ttoto  it  encountered  toco  toot  Jj  »  1,  it  it  moot  toot  ill  devietoric 
stresses  vanish  instantaneously,  o«4  too  volumetric  ttrott  »o  occott  of  T 
vanishes.  Tloit,  too  fiMl  inluif  ttrott  ttoto  tt  Ojj  *  *  1/ J. 

oil  otOor  j  •  0. 

Iwttiflf  tMt  tension  Cutoff  criterion  into  o  vlscoplostlc  fora  tnfort 
the  strottot  or*  released  ot  o  roto  controlled  by  too  fluidity  par*m*ter  y 
rothor  thon  on  Inttontonoout  rolooto.  Accordingly,  It  it  rootonoblo  to 
spoctfy  y  in  the  tension  region  (toy  tj)  ot  o  01 poor  voluo  toon  the  voluc 
of  y  In  the  follure/cop  regions.  moreover,  tlnct  the  tention  cutoff  crl- 
terlon  t roots  volumetric  end  devietoric  ttrott  releotet  independently  the 
viscoplastic  stroln  rote  oust  be  Independently  defined  In  tenet  of  volu- 
metric  and  devlotorlc  strain-rote  components. 

With  the  above  understanding,  the  vlscoplostlc  tention  cutoff  model  it 
defined  by  the  following.  The  plasticity  yield  function  for  tention 
cutoff,  ff.  Is  given  by: 

fT(Jl)  *  Jl  -  T  (20) 
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Hence,  the  “lUtic*  >»eld  surfa <•  (fj  *  Oj  i*  a  stationary  vertical  Ilea. 
J\  *  T.  4t  Shown  in  f  igufe  i.  »**«  t|  0.  tvnyuw*  Cutoff  1%  triggered 
end  the  vlftcopletttc  strain  r«t*  •  *  defined  by 


‘up  *  *T  #^fT}  *  ’G  *(V  *\i  (*l) 

where  fj,  *  /J‘j 

f  J,C  .  I  T 

^  *  *Y  2/j  f  S|1*  2i‘  %}i'  * * n*  ?°?3 

**  j  7 

*.  •  4  >  1 .  1 .  1 .  0,  G.  0 

"  I  (it) 


Equation  21  is  a  Modified  foro  of  the  viscoplesitc  flow  rule 
corresponding  to  Equation  3.  Here  the  first  right-hand-side  tent  contains 
the  vtscoplestlc  strain-rate  for  volumetric  components .  end  the  second  tern 
contains  the  devlatorlc  cowponents.  The  two  tension  fluidity  parameters,  y^ 
and  ^.permit  Independent  control  of  the  volumetric  and  devlatorlc  stress 
release  rates,  respectively. 

As  a  conceptual  illustration,  consider  a  material  that  Is  suddenly 
strained  producing  an  Instantaneous  elastic  stress  state  such  that  fy(Ji)  * 
0.  This  Induces  viscoplastic  flow  (Equation  21)  which  in  turn  releases 
stresses  until  *  0.  When  this  occurs,  we  have  the  steady  state 
condition;  fj  «  fg  *  0,  or  Jj  *  T  and  J'2  *  0,  thereby  satisfying  the  ten¬ 
sion-cutoff  criterion. 

Exact  Solution.  A  deeper  understanding  of  the  tension-cutoff  model  can  be 
achieved  by  obtaining  an  exact  solution  and  studying  Its  behavior.  An 
exact  solution  can  be  obtained  for  a  specified  strain  loading  by  decom¬ 
posing  the  stress  vector  and  elastic  matrix  into  volumetric  and  devlatorlc 
components  and  solving  the  uncoupled  system.  To  this  end,  the  stress  vec¬ 
tor  is  written  as; 
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w 


r 


0  ■  V  ♦  $ 

where •  1  * 
and  i  • 


j  Jj  <  1 »  1 ,  1 ,  0,  0,  0  >  ^  »  volumetric  stress  components 

T 

<  sn,  s22»  S33.  0^2’  °13»  °23  >  *  devlatorlc  stress 


components 


Accordingly,  the  elastic  matrix  is  decomposed  Into  bulk  and  shear  com¬ 
ponents  as: 


K  ♦  6 


(23) 


1 

i 


where,  K  ■  K0 


1110  0  0 
110  0  0 
10  0  0 
sym  0  0  0 
0  0 
0 


bulk  modulus  components 


1 


4  -2  -2  0  0  0 

4-2000 
4  0  0  0 

3  0  0 

sym  3  0 

3 

% 


*  shear  modulus  components 


Now  the  basic  constitutive  relationship,  i.e.,  the  equivalent  of 


Equation  7,  can  be  expressed  in  two  sets  of  equations: 


£  ■  *<£-V 

(24) 

5.  *  -  £,p) 

(25) 

• 

Upon  replacing  cvp  in  the  above  with  tension 

have: 

flow  rule  (Equation  21)  we 

v  =  K  e  -  yt  d>(fT)  K  mT 

(26) 

i 

i 

i 

i 

I  *  G  £  -  Yq  4>(fg)  G  2Jg 

(27) 

i 
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Note,  the  uncoupling  in  the  last  two  equations  is  due  to  the  fact  that  Km,. 
Gi  n»y  «  J). 

If  the  viscous  flow  function  <J>  Is  taken  in  its  simplest  linear  form, 
i.e.  <}>(fy)  *  (Jj  -  T)/fQ,  and  4>(fg)  =  /fg>  then  Equations  26  and  27 
become  the  following  first  order,  linear  differential,  vector  equations: 


v  +  ( 


9K0  yT 


’0 


)  v  =  K  t  +  ( 


3TK0  yT 


)  m7 


sn 

s  +  (  )  s  «  6 

0 


(28) 

(29) 


In  arriving  at  the  above,  use  is  made  of  the  relations,  £  =  K  mT  and 

/p-  9Ko  “  _T 

/J  2  G  mr 

V”"°- 


A  solution  may  be  obtained  for  the  case  of  a  stepped  strain  loading, 

★  ★ 

fc(t)  =  £  h(t),  where  e_  is  any  constant  strain  vector  causing  initial 
elastic  stresses  such  that  Jj  >  T,  and  h(t)  is  a  heavyside  step  function. 
Using  this  loading  in  Equation  28,  we  have  K  e  =  0,  and  the  initial  con¬ 
dition  ^(0)  *  Ke*,  so  that  the  solution  for  volumetric  stresses  is: 

.v(t)  B  (Ke*  *  T/3  my)  exp(-9K0  ^t/fo)  +  T/3  mj  (30) 


Similarly  for  Equation  29,  we  have  G  L  -  0,  and  the  initial  condition  ^(0) 
£  e  *,so  that  the  solution  for  deviatoric  stresses  is: 

i.(t)  =  gc.  exp(-G0  Yg  t/f0)  (31) 


Some  worthwhile  observations  from  these  solutions  are  as  follows: 

(1)  Since  o(t)  *  v^(t)  +i(t),  the  instantaneous  elastic  stress  state 
is  o(0)  =  Kc*+  .Gc*  =  D  e* 

(2)  As  time  increases,  we  eventually  reach  the  steady  state  solution 

v.(®)  =  T/3  mx  and  v(°°)  *  j),  or  =  o22  *  *  T/3,  other  * 

Note,  this  indeed  satisfies  the  tension-cutoff  criterion. 
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(3)  The  exponential  rate  of  stress  release  from  the  Initial  solution 
to  the  steady  solution  is  9K0Yy/f0  for  volumetric  stresses,  and 
GqYg/^o  for  deviatoric  stresses. 

(4)  If  it  is  desired  to  release  volumetric  and  deviatoric  stresses  at 


The  exact  solution  will  be  used  subsequently  to  validate  the  numerical 
solution  algorithm  presented  next. 


Numerical  Solution  Strategy 

In  the  previous  work  effort  (1),  a  numerical  solution  algorithm  was 
presented  for  the  viscoplastic  cap  model  without  tension  cutoff  and  coded 
in  the  program,  VPDRVR.  Here  we  extend  the  algorithm  to  incorporate  the 
tension-cutoff  model.  From  a  programming  viewpoint,  the  incorporation  of 
tension  cutoff  is  straightforward,  only  requiring  modification  to  the 
subroutine  VPLAST  in  the  VPDRVR  program  along  with  the  additional  Input 
data;  T,  Yj  and  Yg  (see  Appendix  A). 

For  completeness,  we  will  review  the  development  of  the  numerical 
algorithm  along  with  a  flow  chart  of  the  complete  model  wherein  the  primary 
objective  is  to  predict  the  stress  reponse  history  from  a  specified  strain 
loading  schedule.  Alternatively,  the  VPDRVR  program  has  the  option  to  pre¬ 
dict  strain  history  from  a  stress  loading  schedule.  However,  the  former 
option  is  directly  suited  for  finite  element  applications. 

Numerical  Approximation.  The  basic  strategy  employs  a  Crank-Nicolson  step- 
by-step  time  Integration  scheme  along  with  a  Newton -Raphson  Iterative  solu¬ 
tion  procedure.  We  begin  with  Equation  7  and  Integrate  over  one  time  step, 
At,  from  time  tn  to  tn+i  to  get  the  incremental  constitutive  relationship: 
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(32) 


Ao  =  D(Ae  -  Ae^) 

where  Ao  =  on+^  -  on  with  on  3  o(t„),  similarly  for  Ac  and  Ac  . 

—  —  —  —  —  n  —  -vp 

All  quantities  at  time  tn  are  presumed  known.  Next,  we  approximate  Ae^  by 
a  one  parameter  Crank-Nicolson  time  integration  scheme  as: 


Ac  =  At  ( (1  -  8  )en  +0  en+1) 
—vp  -vp  -vp 


(33) 


where  0  is  the  adjustable  integration  parameter  in  the  range  0  <_  e  <_  1. 
Choosing  0=0  implies  the  integration  scheme  is  explicit  (simple  forward 

.n 

difference)  so  that  Ae:  is  determined  directly  from  the  known  value  of  e 

-vp  -vp 

at  the  beginning  of  the  time  step.  As  a  consequence,  At  must  be  restricted 
in  size  to  avoid  numerical  instability  (6,7).  Alternatively  choosing  e  > 

0,  the  scheme  is  implicit  since  Ac  is  related  to  the  unknown  value  en+^  at 

-vp  -vp 

the  end  of  the  time  step;  thereby,  requiring  an  iterative  solution  proce¬ 
dure  within  the  time  step.  For  8  >_  0.5,  the  implicit  scheme  is  uncon¬ 
ditionally  stable  so  that  the  choice  of  At  is  governed  by  accuracy,  not 
stability. 


Algorithm  for  Strain  Loading.  Returning  to  the  incremental  constitutive 
relationship  with  Ae^  replaced  by  the  Crank-Nicolson  approximation  and 
using  A a  *  on+^  -  o n  ,  we  rearrange  Equation  32  to  get  the  unknown  quan¬ 


tities  on  the  left  as: 


n"l  n+1  .*  n 

D  a  +  At  0  r 


n+1 


vp 


Ac  -  At ( 1  -  0)^p  +  D'1  a" 


(34) 


Or,  more  compactly,  in  a  symbolic  functional  notation. 
,  i"*1 )  .  a" 


(35) 


where  £  Is  the  vector  of  all  unknown  quantities  at  time  tn+i,  and  is  the 
vector  of  known  quantities  Including  the  specified  strain  increment  Ac  . 
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For  e  >  o.  Equation  34  (or  Eq.  35)  form  a  coupled  set  of  si*  nonlinear 
algebraic  equations  for  the  components  of  '  with  the  understanding  that 

is  to  be  replaced  by  the  appropriate  flow  rule  and  its  associated  yield 
function  depending  on  which  region  of  the  cap  model  is  currently  being 
activated  (i.e.  tension,  failure,  or  cap). 

To  solve  the  above,  a  Newton-Raphson  procedure  is  used  by  expanding  the 
vector  function  £  in  a  limited  Taylor  series  about  a  stress  state  o'  which 
is  some  estimate  of  on+\  and  6o’  is  a  first  order  correction  to  the  esti¬ 
mate,  i.e.  on+1  =  o1  +  601 .  This  leads  to  a  linear  set  of  equations  to 

obtain  the  correction  So1  given  by: 

P*  601  =  -  P1  (36) 


where 

P1  *  D'1  o1 

+  At  6c1 
^vp 

(37) 

and 

ap’ 

P'  =  — — 

—  3  o 

=  o'1  +  Ate 

bo 

(38) 

Here  P’  is  the  Jacobian  matrix  of  the  vector  pi  evaluated  at  a*,  eV. 

—  -vp 

The  iteration  process  is  repeated  with  *  o *  ♦  to  get  a  new 
correction  6o’+^  until  eventually  the  correction  approaches  zero.  Table  1 
summarizes  the  basic  algorithm  (Note  to  start  the  first  iteration  (1*1), 
pi,  P/i,  and  o1  retain  their  values  at  time  tn). 

Upon  reviewing  Table  1  it  Is  evident  that  the  updating  procedure  for 
fi+\  m*^,  e^p1  and  P’^is  dependent  upon  the  current  value  Jj<+^  which 
dictates  what  region  of  the  model  Is  being  activated;  tension,  failure,  or 
cap.  From  a  computational  viewpoint,  the  updating  process  can  be  stream¬ 
lined  by  expressing  all  of  the  plasticity  yield  function,  fc,  fp,  fj, 
and  fg  in  a  general  form  as: 
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T AtU  \  SeUuan  «t*>rUMi 

IkM<X|  MMIM  UUM  •  **#)«§,: 


.  ft  A  •  ft  J*  * 

l.  Given  4t  H«*  1  .  i  •  r  .  |  4ft* 

?.  ring  loop  <i  4  I  (c  n*ii 


(set) 


till  •  oj 


J  ft 


-*9 


i.  Iteration  loop  *  *  1.2.  '*»♦ 

(solve)  P‘*  s^’  *  %*  •  P1 
(update)  £**1  *  J*  • 


s*\ 


,  ft 


i»*l  .■)  »•! 

«,  •  5  - 


i»  r t-t 1 1  .  I  *  f  i«tw>  <«*•«(*  *UJ 


;»r(Jr  :2 )  .  s 


lVJi*  -V  *  **  J» 


(also  cmi^u  ^  »* 


»(*,) 


•  1M 

J’T 

t  * 

Sp 

).  ♦(') 

D»*) 

V- , 

P  * 

p,l4| 

1  £  * 

t '  r 


»r  9 


fV\ 


• 1 


*  I 

tn  c 
Sp 

J^J2 

*2 

4.  Repeat  Iteration  (step  3)  unless  one  of  the  foil 0*1*9  is  satisfied 

(a)  e  »  0,  (explicit  Integration) 

(b)  fn  and  f**1  <  0,  (elastic  dcnaln) 

(c)  |6o^ |  <  tolerance,  (convergence) 

(d)  1  >  imax,  (Iteration  Unit) 

5.  Print  results  and  advance  to  next  load  Increment  (Step  2) 

6.  End 
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*,w: 


where  9iU)  *  -A  *  C  e»p(8l). 


For  the  special  hardening  rule  tor  soils  'not  ro<*sL  i*e  van*  t^r 
steps  are  used  when  loading  the  failure  or  tension  turtles  U*  k). 
except  the  second  step  is  replaced  6* 


-i  +  1 

t  -min 


•  cm  * U .  .... 

!  f*  A  p  y  UX 


where  Xc  -  ♦  R  9  j  C  J  i )  **  1** 

This  completes  tne  1  numeric#!  fo* 


however,  one  last  remark  is  m  0 r*w.  If  *  ♦tlr 


tionship  1$  used,  i.e.#  *  0  Br  wfirrr  3  n  *  f*»*  t  tor*  &*  a*'***. 

*  i  ♦  l 

then  the  update  for  ►.  in  Table  l  must  be  done  it  *ne  1  t  »«• 

i*l  n»l  ^  t*t  i  . ■ '  ♦ 

That  is.  «  t.  -  Vt  where  ^  ’  <r  *  v 


•i*  *  r  ;  +  w  ..  *  r 


is  used  in  the  new  VPDRVR  program  and  is  applicable  f«r  Siot*'  1-ear  #*>£ 
nonlinear  forms  of  the  0  matrix. 


Algorithm  for  Stress  Loading.  For  the  inverse  problem.  -,e..  mpvt  %tr*ss 
loading  with  the  objective  of  determining  the  strain  response,  a 
algorithm  similar  to  that  presented  in  Ta?:e  1  can  b»  eat'ip  established. 

In  point  of  fact  the  stress  loading  algorithm  is  comp ^tat  ♦  oral  ’/  *<utb 
simpler  because  the  associated  Jacobian  matrix  becomes  t*ie  identity  matrix 
as  illustrated  in  the  following  developmer*. 

We  begin  by  rewriting  Equation  34  wttr,  the  unknown  strain  gwartifes  cr 
the  left  and  the  known  quantities  (stress  increment  and  quantities  at  t-me 
tn)  on  right  to  get : 


n+1 
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Atr; 


D  '  .Vt 


.-.til  -  m:;p  * 
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Of,  *4ut«4lt«|t  y.  «»  *  »/WP«ll<  >«Mt  l«Ml 
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n«i 


e*'‘l 

-*P 


It 
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where  the  «0Ct »M  P  and  «J*  are  !►«  left  OM  of  W. 

reaped  »*•!*. 

for  imp  I  »c  it  » ntegra*.  ion  1  *  '  0) .  (  <m<  >«*  V}  >i  a  vtrt  of  tu  nonl  <«Mr 
alg»br#K  equations  tu  b*  suUvJ  if  iur«tt«. 

Af  ternat  laefy  .  *f  »*p  >  >  i  1 1  ‘vf  «%  vied  {»*  0).  fhe  *4«4(  *  wt  are 

t  I f!04f  4*4  do  not  r#m*tr» 

A  pr«t«4k»r*  t«A  It  nUbtiJM  tj  ftptMlAf  M*  | 

limited  'ayfor  series  ah<*»t  (ht  U'H»  State  t*  «*«;*  •  %  M*t  estiamte  of 

n  » 1  f  , 

<  .  and  ’•  M  4  f*r*t  order  .orrection  to  IM  estimate.  1.0. .  i  ♦  ,  •  t>% 

Th*%  leads  to  a  Hdtic  <»l  Of  #01*01*00%  fo*  '<*  ^Irf* 


9 


{*) 


,  —  *  I 

*£ 

Mere  the  Jacobian  natrli  P  becomes  th*  identity  matrla  because  t  CO*  be 

vp 

rtpltctd  by  the  Stress  dependent  flow  rule  so  tf*0t  *  0. 

Based  on  the  abo*e.  it  is  evident  tf»ot  the  algorithm  for  stress  loading 
parallels  the  strain  loading  algorithm  (Table  1)  where  P,  p*,  and  qf  are 
replaced  by  P,  P*.  and  ,  respectively.  Table  3  swart tes  the  stress 
loading  algorithm.  The  procedure  for  updating  the  hardening  paraamters  c  , 

I ,  and  L  is  identical  for  both  algorithms. 

This  completes  the  algorithm  for  stress  loading. 


(44) 

(46) 


when*  9 

9 ' 


Discussion  of  Tension  Cutoff  Algorithm 

In  this  section  we  illustrate  the  performance  of  the  tension  cutoff 
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TABLE  3.  Solution  il^ritm  for  . isoplastic  cap  mode! 
including  tension  culotl  (stress  iMdlnjl 


.  -  n  n  b  n  „«  ,  , 

1.  Given  it  ti«e  t..  *  .  j  ,  t  .  «  .  P  end  tv. 

n  *  r  ~p 

2.  Tine  loop  n  -  l ,  to  hum* 

(set)  q"  1  ft’*  tv  *  •id  -  »)  ♦  tn 

J  Iteration  loop  t  *  I.?.  ima* 


(solve) 

S  c1  * 

at* 

q  -  r 

(update ) 

1  •  1 

t  * 

t  1  *  Sl  ' 

*♦1 

1*1  n-l 

n*l 

L  * 

c  •  0 

> 

1 

< 

1  f  y  ( J|  )  %  *i 

v  T  (also  compute  f^) 

fF(Jr  J,) 

.  l  '  Jj  *  I 

1 

[fc(J,.  Jr 

i)  .  J,  •  L 

m 


•  ♦I 


»♦! 


vp 


.1*1 


•  ijj  (also  compute  if  >  T) 

f»T  «tfT)  mT  ♦  v>  0(fG)  «Iq  .  J,  >  T 


L 


Y>f  ( f  )m  i  Jj  <  T 


L  ‘  V. 

vp 


4.  Repeat  iteration  (step  3)  unless  one  of  the  following  is  satisfied 

(a)  f>  *  0,  (explicit  integration) 

(b)  fn  and  f’*1  <  0,  (elastic  domain) 

(c)  (Sc1 |  <  tolerance,  (convergence) 

(d)  i  >  imax,  (iteration  limit) 

5.  Print  results  and  advance  to  next  load  increment  (step  2). 

6.  End. 
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algorithm  along  with  an  exact  solution.  Also,  a  critique  of  the  model  is 
given  followed  by  some  suggestions  for  future  improvement. 

Illustrative  Example.  Figure  2a  shows  a  tensile  uniaxial  strain  loading 
which  “abruptly1*  jumps  from  0  to  SI  strain  and  held  constant.  Also  given 
in  the  figure  are  the  assumed  elastic  properties  (K0  and  Gq)  and  the  ten. 
sion  cutoff  parameters  >T,  yg,  T  and  f0.  Other  material  properties  asso¬ 
ciated  with  the  failure  and  cap  surfaces  are  immaterial  since  only  the 
tension  surface  is  loaded  in  this  example.  (For  reference,  however,  the  cap 
and  failure  surfaces  were  given  the  properties  for  McCormick  Ranch  Sand 
(1)).  The  viscous  flow  function  Is  chosen  as  $(f)  •  f/f0  so  that  the  pre¬ 
viously  developed  exact  solution  (Equation  30  plus  Equation  31)  can  be  used 
as  a  check. 

In  this  example,  the  exact  solution  for  axial  stress  simplifies  to: 


onU)  -  0.5  exp  (-24t)  +  0.1  (ksi) 

and  lateral  stresses  (o^  =  O33)  are: 


(57) 


o22(t)  =  0.1  exp  (-24t)  +  0.1  (ksi) 


(58) 


These  simple  forms  arise  from  choosing 


deviatoric  stresses  are  released  at  the  same  rate. 


Yj  $0  that  bulk  and 
From  the  above 


equations,  it  is  evident  that  ojj( 0)  *  0.6  ksi  and  022(0)  a  033(0)  ■  0.2 
ksi  so  that  the  instantaneous  value  of  *  1.0  ksi  indeed  triggers 
tension  cutoff  (J^  =  1.0  >  T  *  0.3). 

Figure  2b  shows  the  resulting  stress  histories  as  obtained  from  the 
exact  solution  and  the  numerical  solution  (program  VPORVR).  The  numerical 
solution  overlaps  the  exact  solution  with  less  than  0.2f  error.  Most  of 
this  small  error  is  due  to  the  finite  rise  time  (0.0001  time  units),  used 
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STRESS  REPONSE,  ksi  AXIAL  STRAIN,  €„(%) 


TIME 


Figure  2a.  Uniaxial  strain  loading  for  tension 
cutoff. 


Figure  2b.  Stress  responses  for  tension  cutoff. 
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In  the  numerical  solution,  as  shown  in  Figure  2a,  to  simulate  a  jump 

loading.  It  can  be  observed  that  the  stress  responses  effectively  reach 

steady  state  at  0.2  time  units  after  loading  (on  *  032  ■  033  *  T/3) . 

This  time  period  can  be  increased  or  decreased  by  choosing  yj  higher  or 

9Kq 

lower,  respectively.  It  is  generally  recommended  to  choose  yr  >x — yT  . 

la  laQ  I 

Critique  of  Tension  Cutoff  Model.  Up  till  now,  little  has  been  said  about 
the  rationale  of  the  tension-cutoff  criterion,  l.e.,  is  it  reasonable  to 
assume  that  hydrostatic  tension  (Jj)  by  itself  provides  an  adequate  cri¬ 
terion  for  tension  failure.  From  a  rigorous  viewpoint,  the  answer  is 
generally  no.  However  more  pragmatically,  it  depends  on  the  objective  of 
the  analysis  and,  of  course,  the  type  of  geological  material  we  are  dealing 
with.  Granular  materials,  for  example,  behave  very  erratically  when  one  or 
more  of  the  principal  stresses  are  in  tension.  In  such  cases,  the  ten¬ 
sion  criterion  may  be  as  appropriate  as  any  other  criterion,  particularly 
if  the  analysis  objective  is  to  simply  provide  a  means  of  effectively 
reducing  the  stresses  of  those  elements  experiencing  tension. 

On  the  other  hand,  some  brittle  rock  materials  exhibit  fairly  well 
defined  fracture  planes  when  loaded  in  tension.  Here  more  realistic  tension 
failure  theories  are  available,  such  as,  maximum  principal  stress  theory  or 
Wil 1 iam-Warnke  models  (8).  These  theories  employ  three  Independent  stress 
invariants  (e.g.  Jj,  J ‘ 2  and  J3)  to  describe  the  Initial  tension  cutoff 
surface  and  are  inherently  anisotropic  in  the  post  fracture  analysis.  When 
initial  tension  failure  occurs,  normal  and  shear  stress  components  on  the 
fracture  plane  are  released.  However  on  planes  orthogonal  to  the  fracture 
plane,  stresses  are  still  active.  If  additional  tensile  loading  Is  applied 
such  that  three  fracture  planes  develop  then  all  stresses  are  released  and 
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the  effective  stiffness  in  all  directions  is  zero. 

To  be  sure,  the  tension  cutoff  model  is  not  capable  of  tracking  the 
progressive  fracture  planes  for  brittle  materials.  At  best  It  may  be  said 
that  Jj  model  simulates  a  complete  tension  failure  with  a  small  residual 
hydrostatic  stress  T  which  may  be  specified  as  small  as  desired. 

In  summary  the  Jj  tension  model  may  be  adequate  for  granular  materials 
by  default,  i.e.,  other  tension  cutoff  criterion  have  not  demonstrated  a 
faithful  representation  of  granular  materials  in  tension.  For  brittle 
materials,  the  Jj  tension  cutoff  criterion  is  a  crude  approximation,  and 
more  rigorous  models  are  available.  Nonetheless,  if  the  analysis  objective 
is  to  simulate  loss  of  material  strength  in  localized  areas  of  tension,  the 
model  provides  a  good  engineering  approximation.  Many  soil -structure 
problems,  including  ground-shock  problems,  fall  into  this  category. 

Modifications  for  Viscoplastic  Tension  Damage.  Presuming  that  the  ten¬ 
sion  cutoff  criterion  is  acceptable  for  some  soils  and  rocks,  we  now 
discuss  how  the  viscoplastic  tension  model  could  be  modified  to  represent 
tension  damage  accumulation.  That  is,  limited  experimental  evidence  indi¬ 
cates  that  the  rate  of  tensile  deformation  increases  for  each  loading  cycle 
in  tension  (9).  Conceptually  we  can  conceive  of  this  as  the  progressive 
growth  of  microcracks  which  do  not  heal  during  the  cyclic  loading. 

As  previously  presented,  the  viscoplastic  tension  model  is  insensitive 
to  tension  damage  accumulation  because  regardless  of  how  many  times  it  is 
cycled  into  tension,  the  viscoplastic  strain  rate  remains  proportional 
to  the  fludity  parameter  Yy  and  the  viscous  flow  function  0  (Note  Yq  is 
assumed  to  be  related  to  Yy).  Thus,  if  each  tension  stress  cycle  is  the 
same,  has  the  same  flow  rate  in  each  tension  cycle. 
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To  simulate  the  tension  damage  phenomenon,  two  modifications  are 
suggested;  (1)  a  strain  softening  function  for  the  tension  cutoff  sur¬ 
face,  and  (2)  a  functional  representation  for  Yj.  Both  functions  could 
employ  the  same  history  dependent  measure  for  tension  damage  accumulation. 
As  an  example,  the  strain  softening  function  for  tension  could  be  taken  as: 

Wl>  h)  =  Jl  '  Hh)  (59) 

where  T(£j)  *  T0  exp  (-04  £7)  (60) 

In  the  above  T0  and  aj  are  positive  material  constants  and  £7  is  a  monoto- 
nically  increasing  measure  of  tension  damage  accumulation  (discussed 
subsequently).  Accordingly,  as  tension  damage  accumulates,  fj  Increases 
for  a  given  value  of  thereby  increasing  the  magnitude  of  e^.  Also, 
since  T( c j)  decreases  with  increasing  tension  damage,  the  tension  cutoff 
criterion  is  triggered  at  successively  lower  values  of  in  cyclic 
loading.  This  mimics  non-healing  crack  growth. 

The  functional  modification  for  yj  could  be  taken  as 

Yj(cj)  ~  Y7  exp(ot2  £7)  (61) 

where  yj0  and  <*2  are  positive  material  constants.  Here  YyUyHncreases 
with  tension  damage  thereby  increasing  the  viscoplastic  strain  rate  and  the 
rate  of  stress  release. 

Lastly,  the  measure  of  viscoplastic  tension  damage  accumulation,  £7, 
used  in  both  modifications,  could  be  defined  as  an  accumulation  of  positive 
increments  of  volumetric  viscoplastic  strain  (similar  to  the  cap  hardening 
argument),  i.e.. 
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♦ 


(62) 


eT 


where  e-j.  * 


£VD  ♦  ^VD  +  CVD  *  >  0 

vpll  vp22  vp33  ' 

0,  f  <  0 


(63) 


and  eT  is  an  initial  value  for  tension  damage. 

'o 

The  foregoing  modifications  are  merely  suggestions  to  Indicate  how 
viscoplastic  tension  damage  accumulation  could  be  represented  within  the 
general  framework  of  the  viscoplastic  cap  model.  Although  the  Incorporation 
of  these  modifications  into  a  computational  procedure  is  relatively 
straightforward,  more  experimental  data  is  needed  to  verify  the  validity 
of  these  or  other  possible  forms. 
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PART  II.  PARAMETER  IDENTIFICATION  AND  EXAMPLES 


Introduction 

The  complete  description  of  the  viscoplastic  cap  model  requires  the 
following  identifications;  elastic  parameters  (K  and  6),  failure  surface 
parameters  (A,  B,  and  C),  cap  surface  parameters  (X0  and  R)  along  with  har¬ 
dening  parameters  (W  and  D),  tension  cutoff  parameter  (T),  and  the  viscous 
flow  function  exponent  (N)  along  with  the  compressive  fluidity  parameter  (y) 
and  the  tensile  fluidity  parameters  (Yj  and  Yq).  The  normalizing 
constant  (f0)  should  not  be  viewed  as  an  independent  parameter  and  is 
recommended  to  be  taken  as  f0  =  A. 

For  subsequence  reference,  the  pertinent  functional  forms  employing  the 
above  parameters  are  listed  below  (excluding  tension  cutoff). 


Failure  and  cap  surfaces. 

[fp  =  /J'2  -  (A  -  C  exp(8J1))  ,  T  >  J1  >  L 

fc  =  (J'2  -  (<x  -  L)2  -  (J1  -  L)2)/R2)/fQ  ,  J1  <  L 
Cap  hardening. 


A  =  ln(e/W  +  1 )/D 


£ 


t  . 
e  dt 

0 


min  (dJvp , 

0)  . 

J-j  <  L 

max  (w  . 

0)  . 

J1  >  L 

e  + 

£ 

+  £ 

vpll 

vp22 

vp33 

(soil  only) 


eq  =  W(exp(DXQ)  -  1) 


(64) 

(65) 

(66) 

(67) 

(68) 

(69) 

(70) 
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l  «  X  ♦  R(A  -  C  exp(Bl) ) 


Viscous  flow  rule. 


ivp  *  v*(f)« 


0(0 


(f/fQ)  .  f  >  0 


'0  .  f  <  0 

9'  |  b  ♦  g\,  a 
Constitutive  relationship. 
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3f 

ao 


0  = 


£  - 


in) 

m 

(O) 

(M) 

(75) 

(76) 


All  the  above  form  have  been  developed  and  discussed  In  the  first  part  of 
this  report  as  well  as  in  Reference  (1).  The  tension  cutoff  equations  are 
not  repeated  here  since  the  identification  of  the  associated  parameters  (T, 
Yf  and  yg)  has  already  been  addresed,  i.e.,  T  should  be  selected  as  some 
fraction  of  FCUT  (say  T  =  0.5  FCUT),  and  yg  r  vf/6,  and  >  y  taken  an 
order  of  magnitude  larger  than  y . 

We  now  consider  identi f ication  techniques  for  the  remaining  parameters. 
To  be  sure,  identification  is  a  difficult,  non-unique  process  and  Is  dep¬ 
endent  upon  the  quality  and  quantity  of  experimental  data.  With  regard  to 
identification,  experimental  testing  may  be  grouped  in  two  categories 
“ideal"  experiments  and  "non-ideal"  experiments.  The  former  is  a  complete 
set  of  experiments  explicitly  designed  to  provide  relatively  easy  parameter 
identification  by  controlling  the  loading  schedule  to  permit  separate  iden¬ 
tification  of  viscous  and  plastic  (steady-state)  responses.  Conversely, 
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the  "noo- !  d*a 1  “  *xp*r l»*nt S  *r*  rlr  *nj  v*  Jk,  ♦"•cl  r*»,»><Ae  aiint 

information  on  uncoupling  viscous  and  ii'astK  >UI* 

For  the  non-ideal  Uitf,  parameter  id*nt i » ical ion  is  rr«&aD‘>  5til 
by  a  trt*l  and  error  approach  usinj  th»  program. 

Parameter  idem I f  icat loo  for  Doth  in*  i<s*al  and  non- ideal  <ase>  H 
discussed  in  subsequent  sections.  unfortunately  *ci!  of  *.r«  «un»j 
experimental  data  for  toil  and  r©«,»s  fails  »mo  in*  'non-ideal*  taiepory. 
To  circumvent  this  problem.  we  alll  dafin*  an  ideal  set  of  experiments 
along  with  a  hypothetical  set  of  results  to  illustrate  the  'ideal*  iden¬ 
tification  procedures.  Hopefully  this  will  serve  as  a  template  for  future 
experimental  work.  For  the  non-ideal  cases,  existing  experimental  data 
will  be  used  to  illustrate  parameter  identification  by  a  trial  and  error 
process  using  the  VPDRVR  program. 

Before  these  identification  procedures  are  presented,  it  is  well  to 
discuss  some  behavioral  aspects  of  the  vtscoplasHc  cap  model  along  with 
the  influence  that  various  parameters  have  on  the  model’s  response. 

Model  Behavior  and  Parameter  Influence 

The  behavior  of  the  viscoplastic  cap  model  is  best  visual t/ed  by  con¬ 
sidering  the  plasticity  surfaces  In  *P«ce  as  shown  m  Hgjr*  >a 

Here  the  Initial  cap  setting,  X0j  defines  the  Initial  elastic  region  shown 
by  the  shaded  area.  To  begin  with,  we  suppose  a  compressive  stress  state 
"oa"  is  suddenly  applied  and  held  constant  such  that  fc(’^,  XQ)  9  0  but 

within  the  failure  surface  as  illustrated  in  Fig.  3a.  Viscoplastic  flow 
will  commence  with  a  rate  proportional  to  >$(fc)m.  fq.  72.  At  the  same 
time,  the  cap  location  X  will  begin  to  move  from  Its  Initial  location  X0 
toward  the  location  Xa  so  that  as  time  passes  fc(oA>  x)  (end  hence  evp) 
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Figure  3b.  Influence  of  parameter  R  on  viscoplastic  st*a'«s. 
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to  both  model  shapes.  As  before  if  a  sires'.  stale  is  s-dOeelj  imposed, 
both  cap  surfaces  wilt  move  from  tneir  initial  locations  «ro  eteniwaiij 
reach  a  steady-state  position  containing  me  stress-state  *  as  show*. 
During  this  movement,  the  relative  distribution  ot  •olumetric  ana 
deviator'c  viscoplastic  st'atrs  is  -,jvwiirJ  t*e  u.t»arj  n»f#a ‘  at 
Consequently  the  vertical  ell’pse  (w  1 >  d'stnbutes  a  grea ter  j  roper  turn 
of  its  scoplast  ic  strain  tt  tne  volumetric  components  ’.n4n  does  me  h©ri> 
ZOntdl  ellipse  (ft  •  11.  Thus,  thi-  major  purges**  of  tne  ;a'amrter  *  is  to 
control  the  relative  distft  utior:  ’  .ui'ivt''  t*.  ar  .  *t*.  at  <  '.  *•  •.  *v 

strain  components. 

Cap  hardening,  which  locates  tne  current  ^c%’'ton  ot  1.  **•%>  o»s  tne 
parameters  W  and  D  along  with  the  accumul  at  wiscopi«st  *c  1U1  f  harden1*' 
measure  as  given  by  tquation  66.  However  tne  mfl.enoe  >f  a  and  t  is 
much  more  complicated  than  this  equation  initially  suggests  because  o'H«M- 
tely  is  also  dependent  on  W  and  l!  through  the  v'SvOp'ast  it  t‘«ow  -/e. 
Nonetheless,  if  Equation  bo  is  plotted  m  norqireng * ona I  * ora  as  shown  ir, 
Figure  4,  it  can  be  used  to  great  advantage  in  understand* r ,  tne  hardenmj 
behavior,  as  well  as,  helping  to  quantify  the  W  and  0  parameters. 

Listed  below  are  some  useful  insignts  with  regard  to  the  harden *rq 
relationships  in  Figure  4. 

(1)  D  (units  of  inverse  stress)  and  U  (units  of  strain)  are  positive 
constants  so  that  DX  and  -/W  are  dimensionless  hardening  coor¬ 
dinates.  At  any  instant  in  time  DX  and  > /M  are  related  by  the 
graph  in  Figure  4. 

(2)  Assuming  the  initial  value  of  X  =  X0  is  negative  (which  also  gives 
fo<0)  then  all  subsequent  values  of  X  and  c  are  negative  such 
that  X0  ^  -"and  tQ  _>  t  _>  -  W.  Thus  W  represents  the  maximum 
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absolute  value  of  volumetric  vjscoplastic  strain  accumulation  that 
can  ever  be  achieved  under  any  loading. 

(3)  When  l/W  *-1.0  or,  equivalently,  when  ;DX;  >  3.0  (approximately) 
we  see  that  a  small  change  in  At  -  eAt  produces  a  large  change  AX  « 
XAt.  As  a  consequence  X  moves  its  position  (hardens)  very  rapidly 
toward  the  steady-state  position  so  that  the  net  response  is 
almost  entirely  elastic  (i.e.,  there  is  no  time  for  viscoplastic 
strains  to  accumulate).  This  is  why  the  cap  model  produces  an 
asymtot  ical  ly  increasing  stress-strain  curve  in  hydrostatic 
compression,  in  effect,  a  pure  elastic  slope  in  the  limit. 

Based  on  the  above  insights,  experience  has  shown  that  it  is  usually  best 
to  quantify  the  parameter  D  in  conjunction  with  X0.  That  is,  X0  is  first 
chosen  to  establish  the  size  of  the  initial  elastic  domain,  then  0  is  cho¬ 
sen  such  that  |DX0|  •-  3.0  in  order  that  plastic  hardening  will  be  effec¬ 
tive.  Typically  |DX0'  <  0.5  has  been  used  in  this  study.  Having  made  the 
selection  for  D,  the  parameter  W  may  be  used  for  control.  Increasing  W 
results  in  increased  amounts  of  viscoplastic  strain,  and  conversely, 
decreasing  W  decreases  the  amount  of  viscoplastic  straining. 

Identification  Guidelines 

Summarized  below  are  the  behavioral  aspects  and  parameter  iden¬ 
tification  guidelines  for  the  viscoplastic  cap  model. 

Elastic  Parameters.  The  bulk  modulus  K  and  shear  modulus  G  control  the 
elastic  response.  These  parameters  are  best  determined  from  unloading 
tests  and  are  restricted  such  that  K  >  2/3G.  Nonlinear  forms,  K  =  K ( J j )  and 
G  =  G(J2)»  may  be  used  if  warranted. 

Failure  Surface.  The  failure  surface,  denoted  by  fp  in  Equation  64,  is 
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defined  by  parameters  A,  B,  and  C  and  forms  a  static  yield  surface  fp  =  0 
along  which  the  cap  surface  can  move.  If  a  constant  stress  state  Is 
imposed  such  that  fp  >  0,  steady-state  failure  conditions  will  result 
wherein  the  deviatoric  components  of  tvp  will  continue  to  increase  at  a 
constant  rate.  The  parameter  A  is  the  maximum  height  of  the  static  failure 
surface  and  A-C  is  the  height  at  J]  =  0,  hence  A  >  C.  The  parameter  B 
controls  the  curvature  of  the  static  yield  surface,  as  B  incrases,  the  rate 
of  curvature  increases. 

Cap  Surface.  The  cap  surface,  denoted  by  fc  in  Equation  65,  is  an  ellipse 
quadrant  whose  shape  is  governed  by  the  parameter  R  and  whose  initial  loca¬ 
tion  is  set  by  X0.  Choosing  R  >  1  implies  a  horizontal  ellipse  whereas  R  < 

1  implies  a  vertical  ellipse.  By  increasing  R  the  outward  normal  of  the 
cap  surface  leans  more  toward  the  deviatoric  direction  (  /j*2  )  so  that  a 
greater  proportion  of  evp  is  weighted  toward  the  deviatoric  components  as 
opposed  to  the  volumetric  components.  In  a  typical  triaxial  stress  loading, 
for  example,  the  effect  of  increasing  R  is  to  increase  the  total  axial 
strain  tjj.  As  a  first  guess,  use  R  =  1.0  for  trial  and  error 
identification. 

Setting  X0  to  a  large  negative  value  provides  a  large  initial  elastic 
space  which  is  often  a  suitable  representation  for  pre-consol i dated  soils 
or  rocks.  However  remolded  soil  specimens  are  usually  best  represented 
with  an  initially  small  elastic  domain  in  which  case  X0  is  assigned  a  small 
negative  value  near  the  origin.  Once  X0  is  selected,  e0  is  given  by 
Equation  70,  and  L0  is  given  by  Equation  71. 

Cap  Hardening.  The  parameters  D  and  W,  Equation  66,  control  cap  hardening 
along  with  the  accumulated  volumetric  viscoplastic  strain  measure  e  , 
Equations  67-70.  Experience  has  shown  it  is  best  to  set  0  <  1/2  |  X0|  "1 
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and  use  W  to  control  the  cap  hardening.  To  wit,  increasing  W  retards  the 
cap  movement  thereby  increasing  the  accumulation  of  evp.  As  |e/Wj 
approaches  1.0,  the  rate  of  cap  movement  X  becomes  infinite,  and  hence,  the 
response  becomes  more  elastic. 

Viscous  Parameters.  The  cap/failure  fluidity  parameter  y  linearly  controls 

the  viscoplastic  strain  rate  e  vp  via  Equation  72.  The  exponent  N  also 

influences  the  rate  but  in  a  nonlinear  fashion  dependent  on  the  current 

N 

value  of  f /f o »  ^vp  =  V-J>(f/fo)  ?!•  Typically  N  =  1  is  used,  and  £vp  is 

controlled  by  y  unless  sufficient  experimental  data  is  available  to  quan¬ 
tify  both  parameters.  The  tension  fluidity  parameters  are  usually  chosen 
on  an  arbitrary  basis  since  experimental  data  in  tension  is  usually 
lacking.  It  is  recommended  to  choose  yj  an  order  of  magnitude  larger  than  y 
and  set  Yq  =  9  K/G  Yy. 

Tension  Cutoff.  The  tension  cutoff  parameter  T  is  also  usually  chosen  on 
an  arbitrary  basis.  However,  it  is  limited  to  a  small  practical  range,  0 
T  FCUT,  where  FCUT  is  the  intersection  of  the  failure  surface  with  the  Jj 
axis.  Typically,  we  set  T  =  FCUT/2,  but  we  must  check  that  T  >_  L0  if  L0  is 
positive. 

The  foregoing  guidelines  are  used  extensively  in  the  following  iden¬ 
tification  procedures. 

Parameter  Identification  Using  Ideal  Experimental  Data 

In  the  following,  an  ideal  experiment  is  defined  and  hypothetical 
results  are  constructed  to  illustrate  an  identification  procedure. 

Loading  Schedules.  A  conventional  triaxial  testing  apparatus  is  used  to 
conduct  the  ideal  test.  We  require  at  least  one  purely  hydrostatic  test 
and  a  minimum  of  three  triaxial  tests  with  successively  increased  confining 
pressures.  Figure  5a  illustrates  the  stress  paths  of  these  tests  in  terms 
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of  oi  and  03  where  H  is  the  hydrostatic  test  (01=  03)  and  Tj,  Tg,  T3  are 
the  triaxial  tests.  Each  test  employs  a  virgin  material  specimen,  and  as 
illustrated,  we  require  an  unload-reload  cycle  on  the  hydrostatic  test  and 
one  of  the  triaxial  tests. 

The  time  loading  schedule  for  the  hydrostatic  tests  is  illustrated  in 
Figure  5b.  Here  a  sequence  of  pressure  increments  are  applied  such  that 
each  increment  is  maintained  until  steady-state  conditions  are  observed 
(i.e.  no  volume  change,  w  =  0).  Similarly,  the  time  loading  schedule  for 
each  triaxial  test  is  a  sequence  of  axial  stress  increments  as  shown  in 
Figure  5c.  To  start,  we  require  the  initial  hydrostatic  loading  to  be  at  a 
steady-state  condition  before  the  first  axial  load  increment  is  applied. 

Each  axial  load  increment  is  applied  and  held  constant  until  a  steady-state 
condition  cj  =  0  is  observed  for  axial  strain  after  which  the  next  load 
increment  is  applied.  Eventually  a  failure  condition  is  observed  ej  ■ 
constant  (or  increasing  in  rate)  which  terminates  the  test.  In  both  types 
of  tests  an  unload-reload  cycle  is  conducted  from  a  steady-state  position. 
Response  Data.  For  the  hydrostatic  test  we  require  steady-state  data  of 
pressure  vs.  volumetric  change,  as  well  as,  a  time  history  plot  of  volu¬ 
metric  change  vs.  time.  For  each  triaxial  test  we  require  steady-state  data 
of  shear  stress  -  03  vs.  axial  and  lateral  strains.  As  the  first  step 
in  parameter  identification,  we  will  use  steady-state  data  to  identify  the 
plasticity  parameters.  To  this  end,  the  hypothetical  stress-strain  responses 
are  plotted  in  Figures  6a  and  6b. 

Identification  Procedure.  The  basic  strategy  for  parameter  Identification 
is  to  determine  elastic  parameters  from  unloading  data,  plasticity  parameters 
from  steady-state  data,  and  viscous  parameters  from  time  history  response 
data.  This  Is  illustrated  in  the  following  steps. 
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PRESSURE  -J|/3 


Figure  6a.  Pressure  vs.  volume  change  at  steady  states. 
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(1)  Elastic  parameters  (K,  G).  The  bulk  modulus  K  is  the  slope  of  the 


hydrostatic  unload-load  curve  (Figure  6a)  and  Young's  modulus  is 
the  slope  of  the  triaxial  unload-load  curve  (Figure  6b).  Thus,  G 
=  3KE/(9K-E). 

(2)  Initial  elastic  domain  (X0).  The  initial  elastic  domain  as 
controlled  by  X0  can  be  determined  by  transferring  the  bulk  modu¬ 
lus  slope  K  to  the  origin  of  Figure  6a  and  finding  the  presure  P0 
at  which  the  slope  departs  from  the  data  curve.  Thus,  X0  =  3P0. 

(3)  Failure  surface  (A,B,C).  The  failure  points  for  tests  Tj,  T?,  and 
T3  (Figure  6b)  represent  the  maximum  steady-state  shear  stress 
obtainable  for  each  confining  pressure.  These  three  points  may  be 
plotted  in  Jj,  / J 1 2  space  with  Jj  =  oj  +  2^3  and  Sj'2  *ioj  -  03I//3 
as  shown  in  Figure  7.  Accordingly,  the  parameters  A,  B,  and  C  may 

be  determined  by  the  failure  condition  fp  =  0,  (i.e.,  /j‘2  =  A  -  C  exp(BJj) 
for  the  three  data  points.  This  may  be  done  graphically  by 
establishing  the  values  for  A  and  C  as  shown  in  Figure  7  and 
using  the  above  equation  to  compute  B.  Alternatively,  a  least- 
squares  error  fit  could  be  used  to  get  A,  8,  and  C  simultaneously. 

(4)  Cap  hardening  (W,  0).  Since  hardening  is  controlled  by  the  volu¬ 
metric  viscoplastic  srain  measure  £,  the  steady  state  hydrostatic 
data  (Figure  6a)  is  sufficient  to  determine  the  hardening  parame¬ 
ters  W  and  D  by  using  Equation  66.  To  this  end,  we  note  that  the 
steady-state  hydrostatic  data  implies  X  =  Jj  (i.e.,  X  coincides 
with  Jj  at  steady-state)  and  e  =  eQ  +  wvp  (i.e.,  wvp  is  always 
compressive  so  that  all  increments  add  to  e).  Since  wvp  =  w  -  J ^ /3k 
we  can  say  c  =  e0  +  w  -  Jj/3K,  where  e0  is  defined  by  Equation  70. 
Substituting  the  above  relationships  for  X,  c  ande  0  into  Equation 
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Figure  7.  Failure  surface  plot. 
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W 


66  we  have 


W(exp(DJi)  -  exp(DX0))  -  w  -  Jj/3K 

In  the  above  K  and  X0  are  known  constants  (Steps  1  and  2}  and  w 
and  Jj  are  steady-state  data  points  in  Figure  6a.  Accordingly, 
the  unknown  parameters  W  and  D  are  to  be  determined  to  best 
satisfy  the  above  equation  for  all  data  points.  The  easiest  way 
to  do  this  is  to  choose  I)  ^  1/2  i  X0i  ^  and  directly  solve  for  W 
at  several  (Jj,  w)  data  points.  If  each  of  the  W’s  so  determined  is 
not  approximately  the  same,  make  a  small  adjustment  in  D  and  try 
again. 

(5)  Cap  shape  (R).  The  strategy  to  determine  R  is  illustrated  in 

Figure  8  where  a  particular  failure  data  point  (T2)  is  chosen  with 
known  coordinates  Jj*,  O?*.  At  this  point,  we  have  L  =  Jj*,  and 
our  objective  is  to  determine  the  location  X  so  that  R  is  given  by 
R  =  (L  -  X)/  >/j2*.  To  get  X  we  use  the  hardening  function.  Equation 
66,  in  which  W  and  D  have  already  been  determined,  so  our  problem 
is  to  find  £  at  the  failure  point  Jj*,  >^2*.  This  may  be 
achieved  by  adding  the  volumetric  viscoplastic  strain  from  the 
hydrostatic  test  evaluated  at  with  the  additional  volumetric 
viscoplastic  strains  from  triaxial  test  T2  (Figure  6b).  Thus  the 
computational  steps  are  as  follows: 
w  =  w(J]*)  +  (Aej  +  2Ac3) 
wVp  =  w  -  Jj*/3K 

£  =  £q  +  Wyp 

X  =  ln(e/W  +  1)/D 


(volumetric  strain) 
(viscoplastic  volumetric 
strain) 

(hardening  strain) 

(cap  X  associated  with 
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lb) 


K  -  (Jj*  -  *)/  »v  /  of  * 

This  process  may  br  repeat**!  fur  ^  ai^  f«'!w*r  ,  c '■  M  s 
js  desired  and  an  a^eragr  ‘  <*r  ?*f  *  .  «r  !x  Jk-t  crwinre. 

Vi  scous  Parameter  s  • ,S; .  . . 1  -  ,  •  ^  Vj  r^iu  cut  #  *^r  ,  niMt 

ters  ’  and  N  can  :>**  detwm  t  n#-  ?  f«  ••  ■  •  «*  •  >d' v!  J*.  *  t  fri*  t » 
ma  t  c  h  1  r  9  v  1  sc  op  1  a  s  t  '  ^.irnr  •  j  t  «*  »  » i  !  t  t  hr 

viscop!artK  f )  ow  rule,  ^4!  io~  ' .  .  !*,  br;**  wit*.  *r  j.  1  <M  w4 

vs.  time  as  illustrated  •  r  9  »hrr*  thr  1 1**  s?*ti^s  !j* 

t£,  etc.  correspond  tr  the  hydrostatic  ioadtrK  %.hed^  e  in  Mgwfr 
5b.  At  the  beginning  of  each  creep  pha»**.  e.  .  time  -  t  ^  'f' 


Figure  9,  the  volumetric  viscoplastu  strain  ra!«*  .^p*  .  .  *r  !*• 

measured  as  the  tangent  to  th«»  1 1  nr  htstutj  p’ot.  Vuc^  rat* 

values  wvp  ( t  i )  ma>  *>e  obtained  for  each  time  station  \  y  *  t «  ,  tj# 

tj,  ...  terminating  a  load  increase.  Equating  t**is  data  tc 

the  volumetric  viscoplastic  flow  rule  we  have 

•  * 

w  (t . ;  f  <  f  f  ’k<v > ,  1 

vp'  1  c  0  v  I 

where  mj  ?(Jj  -  l)/f0  *<' 

It  is  understood  that  f-  and  mj  are  to  be  evaluated  a*  t'me  t . 
Evaluating  fc  and  mj  is  laborious  but  straight  forward  as  follows, 
first  compute  (tj)  -  •  0  ♦  wvp(tj),  second  compute  X*t.r  frc*r 
Equation  66,  and  third  compute  l(t*)  from  equation  7|.  Havmc  I 
and  L  along  with  Jj(tj),  fc  and  mi  may  be  directly  evaluated  s^nce 
all  plasticity  parameters  are  known  and  J'^  2  0. 

Now  for  each  wVp(t<j)  iota  point  everything  is  known  in  thf 
above  expression  except  y  and  N.  These  parameters  can  be  deter¬ 
mined  by  a  least-square  error  technique.  Or,  more  simply,  choose 
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VOLUMETRIC  VISCOPLASTIC  STRAIN 
Wvo  •  W-  J/SK 


W 


* 


TIME 


Mcjure  9  vtscop* ast*c  *%*■*»*  ;^r\ 


N  =  1  ana  determine  vabje-  u  It  »r  i  '-u-VuleO  *  ‘s  ar* 

not  in  good  agreement  iJ‘u^  N  a"  :  t'y  d^'n. 

This  completes  the  ideal  lCent  •  f  Uat  tor*  p'v^ess,  we  ?«*•*!  consiaer  t  r  5 
and  error  methods  with  the  VPOK^h  program. 

Parameter  Ident  1  f  icat  ion  Usi  nq_  Non  -  Idea  ’  ^  .  Data 

Three  distinct  sets  ot  e»p«'ri  :.rnt«*l  ale  co'*s  tdenrj  #  /  ift*  ^uf*.-usr  uf 
identifying  the  vis^opiastic  parameters  by  a  t  m  a  ’  and  e'fur  procedure 
using  the  V  PDKVR  program.  The  three  experiments  represent  *  ^arg»*  of 
geologi  cal  materia  : nard  'imestone.  soft  sednent  sry  rou,  and  we  i  1  -gr<*v1»<3 
sand.  Further,  the  manner  of  loading  and  loading  rates  are  si  *r  •  *  team '/ 
different  between  each  experiment,  thus,  tnis  study  not  on  ^  illustrates 
tne  parameter  identification  process,  hut  also,  d^mo/ntra t^s  the  rapab*.’- 
ties  ana  limitations  of  tne  viscoplastic  cap  mode’,  fcacn  experirr^^t  arj 
corresponding  parameter  fit  is  discussed  in  turn. 

Limestone  in  Triaxial  Stress.  A  rather  elaborate,  nonstandard,  t nana1 
test  experiment  on  specimens  of  Solenhofen  Limestone  was  conducted  by 
Robertson  (10)  to  measure  the  axial  strain  history  *  jj  resulting  from,  a 
variable  axial  stress  loading  sequence.  Details  of  tne  testing  apparatus 
and  experimental  program  are  somewhat  involved  and  are  not  repeated  here. 
Instead,  we  simply  identify  the  stress  loading  history  (Figure  10)  for 
Robertson's  specimen  number  S-90  whicj^i  is  considered  in  this  study.  As 
shown  in  Figure  10,  an  initial  triaxial  stress  state  is  rapidly  imposed  (*jj 
=  96.1  k  s i ,  0  2  2  s  J 33  “  44,1  k$i).  Thereafter,  the  lateral  stresss  are 
maintained  constant,  and  the  axial  stress  is  intermittently  step  loaded  at 
time  =  7.2,  12.9,  and  22.8  Mloseconds.  After  each  step  loading  including 
the  initial  loading,  Ujj  decreases  by  some  amount  due  to  the  nature  of  the 
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Cap  Hardening: 


Bulk  •  3,300  ksi  ;  shear,  6(J2')  * 

91  ( I  ♦  21  exp  (-0.0012  J2'  ))  ksi 

V  ^  J|  *  )  *  ^2  “  ( 1.0  -  0.25  J,)  ksi 

R  •  2.4 ;  fQ  -  1.0  ksi  ;  XQ  •  -212.0  ksi 
W  •  0.55  ;  00  -  0.0024  ksi”1 


Viscous  Flow  Function: 

N  ■  I 

r  ■ 


»  f0 
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1.0  ksi 
10* 
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-i 


Figure  10.  Triaxial  stress  loading  schedule  and  model 
parameters  for  limestone. 
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hydraulic  testing  apparatus.  Although  the  magnitude  of  these  decreases 
were  reported,  their  time  history  was  not.  Accordingly,  the  linearly 
decreasing  functions  following  each  jump  in  Figure  10  are  approximations. 

Axial  strain  measurements  were  recorded  before  and  after  each  loading 
step,  providing  a  data  base  for  attempting  to  "curve  fit"  the  viscoplastic 
cap  model.  Since  the  experiment  represents  a  consecutive  sequence  of 
loadings,  "curve  fitting",  in  this  case,  is  quite  difficult  because  the 
accumulated  strain  depends  upon  the  entire  loading  history  and  the  strain 
hardening  parameter  c  controlling  the  cap  movement. 

Figure  11  shows  strain  history  data  points  along  with  a  viscoplastic 
cap  model  representation  producing  a  fairly  good  correlation.  Since  the 
viscopiastic  model  was  driven  by  the  triaxia!  stress  loading  schedule  in 
Figure  10,  the  stress  loading  algorithm  was  used  in  the  VPDRVR  program.  The 
final  parameters  for  the  viscopiastic  cap  model  are  also  shown  in  Figure  10 
and  were  largely  determined  by  trial  and  error,  discussed  next. 

Isotropic  elastic  parameters,  bulk  modulus  and  shear  modulus,  were 
determined  by  best  fitting  the  instantaneous  jump  responses,  i.e.,  no 
viscopiastic  flow  was  assumed  to  occur  during  the  loading  jumps.  This  was 
best  matched  by  a  constant  bulk  modulus  and  a  variable  shear  modulus  mono- 
tonically  decreasing  with  J ‘ 2  (Fig.  10).  The  lailure  surface  was 
simplified  to  a  standard  Drucker-Prager  form  and  the  initial  cap  surface, 
shaped  as  a  horizontal  ellipse  R  =  2.4,  was  located  well  into  the 
compression  region  by  setting  Xq  =  -212.0  ksi.  The  motivation  for  this 
initial  setting  was  to  provide  a  large  elastic  region  so  that  the  initial 
jump  loading  did  not  cause  excessive  viscoplatic  flow  in  accordance  with 
observations.  Also  it  insured  the  viscopiastic  flow  in  accordance  with 
the  cap  surface  (J^  <_  L)  throughout  the  loading  schedule. 
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AXIAL  STRAIN,  -€„(%) 


Figure  11.  Axial  strain  response  and  viscoplastic  cap  model 
representation  for  limestone. 
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In  accordance  with  previously  presented  guidelines,  the  hardening 
parameter  D  was  set  at  1/£|X0|’^.  The  hardening  parameter  W  as  well  as  the 
fluidity  parameter  Y  were  adjusted  by  numerical  experimentation  to  best 
match  the  data.  No  attempt  was  made  to  vary  the  exponent  of  the  viscous 
flow  function  and  was  set  at  N  =  1. 

Although  it  is  not  claimed  the  model  parameters  chosen  here  are  repre¬ 
sentative  of  the  limestone  material  in  any  loading  environment,  we  do 
assert  that  the  representation  in  the  range  considered  is  fairly  good.  It 
must  be  said,  however,  that  Robertson's  data  extended  beyond  the  range 
presented  here,  and  it  was  observed  that,  as  the  axial  load  increased,  axial 
strains  were  increasing  at  an  ever  Increasing  rate.  Such  behavior  may  be 
attributed  to  strain  softening  which  is  not  within  the  capabilities  of  the 
current  viscoplastic  model.  This  is  illustrated  in  another  manner  in  the 
next  experiment. 

Sedimentary  Rock  in  Triaxial  Stress.  The  viscoplastic  yielding  of  soft 
sedimentary  rock  samples  was  Investigated  by  Akai,  et  al.  (4).  Their 
experiments  consisted  of  standard  triaxial  tests  on  cylindrical  samples  of 
a  porous  tuft  described  as  an  ideal  soft  sedimentary  rock. 

The  data  considered  here  is  for  four  separate  creep  tests  all  with  the 
same  confining  pressure  and  different  axial  loads.  Each  axial  load  is 
rapidly  applied  and  held  constant  for  the  duration  of  the  creep  test,  up  to 
8,000  minutes.  Figure  12  defines  the  imposed  stress  states  for  each  of  the 
four  tests  along  with  the  initial  cap  model  setting  and  the  model  parame¬ 
ters  used  for  this  study.  The  measured  strain  history  data  (deviatoric 
strain,  Cjj)  reported  by  Akai  is  shown  by  data  points  in  Figure  13  along 
with  the  viscoplastic  model  representations  shown  with  solid  lines.  Here  it 
is  observed  that  reasonable  correlation  with  the  data  was  achieved  in  the 
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Failure  Surface 


State 


°axial 
( ps  i ) 

lateral 
( psi ) 

497.9 

71.13 

668.61 

71.13 

711.29 

71.13 

739.74 

71.13 

Initial  Cap 


Elastic  Moduli : 

K  *  125,000  psi  ,  G  -  60,000  psi 

Failure  Surface: 

fF  *  ^7  ~  <  275.0  -  0.863  J,  ) ,  p«i 

Cap  Surface: 

8  *  0.35 ,  XQ*  —800  psi 

Cop  Hardening: 

0  *  0.00078  (psi )"',  W*  0.25 

Viscous  F low : 

X*  0.5  x  KT®  { min)"*  ,  N  -  1.6, 

f0  »  275  psi 

Figure  12.  Stress  states,  initial  cap  setting  and  parameters 
for  soft  sedimentary  rock. 

primary  and  secondary  creep  range,  but  not  in  the  tertiary  creep  range. 

This  will  be  elaborated  further  after  a  brief  discussion  on  the  parameter 
identification  procedure. 

The  elastic  properties  were  determined  by  assuming  no  viscoplastic 
flow  occured  while  each  axial  load  was  imposed  so  that  the  initial  strains 
were  elastic.  As  shown  in  Figure  12,  the  initial  cap  setting  was  taken 
well  into  the  compression  range  with  X0  =  -800  psi  and  R  =  0.35  along  with 
a  standard  linear  Drucker-Prager  failure  surface.  The  motivation  for  these 
choices  were  due  to  the  observation  that  the  strain  response  data  exhibited 
continued  elastic  behavior  for  test  1,  creep  and  then  steady-state  behavior 
for  test  2,  and  creep  and  then  steady-state  "failure"  followed  by  tertiary 
creep  for  tests  3  and  4.  Accordingly,  the  initial  cap  setting  was  located 
between  stress  states  1  and  2  to  insure  an  elastic  response  for  test  1. 

The  failure  surface  was  located  slightly  above  stress-state  2  to  achieve 
steady-state  response  for  test  2  and  below  stress-states  3  and  4  to  achieve 
steady-state  failure.  Of  course  there  is  nothing  unique  about  the  par¬ 
ticular  parametric  values  chosen  to  accomplish  this  initial  setting.  In 
accordance  with  previous  guidelines,  the  hardening  parameter  0  was  taken  as 
a  fraction  of  |X0|~1  and  the  remaining  parameters  W,  Y  ,  and  N,  shown  in 
Figure  12,  were  chosen  by  numerical  experimentation  with  VPDRVR  program. 

Returning  to  the  model's  performance  shown  in  Figure  13,  we  observe 
the  elatic  response,  test  1,  and  the  steady-state  viscoplastic  response, 
test  2,  are  well  correlated  with  the  experimental  data.  In  test  3,  the 
model  correlates  fairly  well  with  primary  and  secondary  creep  data. 

Primary  creep  is  the  early  portion  of  the  curve  with  decreasing  strain 
rates  (cap  movement)  and  secondary  creep  is  a  constant  strain  rate 
(steady-state  failure).  The  last  two  data  points  in  test  3  exhibit  ter- 
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tiary  creep,  I.e.,  increasing  strain  rate,  which  is  not  represented  by  the 
viscoplatic  cap  model.  Test  4  also  exhibits  tertiary  creep  beginning 
almost  immediately  after  the  primary  creep  phase.  Again,  tertiary  creep  is 
not  represented  by  the  viscoplastic  cap  model. 

From  these  comparisons  we  conclude  that  the  present  viscoplastic  cap 
model  is  capable  of  simulating  elastic,  primary  creep  and  secondary  creep 
behavior  but  not  tertiary  creep.  The  inability  of  the  viscoplastic  model 
to  simulate  tertiary  creep  is  not  a  question  of  readjusting  the  parameters, 
but  rather,  it  is  an  inherent  limitation  of  the  functional  forms  defining 
the  model,  i.e.,  the  present  model  can  only  respond  with  constant  strain 
rates  once  the  steady-state  failure  condition  is  reached  (e.g.,  recall  Fig. 
3a).  One  way  of  overcoming  this  limitation  is  to  introduce  strain  soften¬ 
ing  into  the  hardening  function  such  that  after  e  has  grown  (hardened) 
to  a  specified  level,  a  softening  function  is  activated  shrinking  c,  and 
hence  increasing  the  strain  rate  as  the  cap  retracts.  Another  approach 
would  be  to  redefine  the  fluidity  parameter  in  a  functional  form  dependent 
on  r  .  This  idea  was  discussed  at  the  end  of  Part  1  of  this  report. 

Sand  in  Uniaxial  Strain  with  Variable  Load  Rates.  We  now  consider  the 
last,  and  perhaps,  the  most  significant  experimental  test  for  evaluating 
the  performance  of  the  viscoplastic  model,  as  well  as,  identifying  the 
model's  parameters.  This  rather  ingenious  experimental  test,  conducted  at 
the  Army's  Waterways  Experimental  Station  (11),  was  undertaken  to  directly 
assess  the  effect  of  loading  rate  on  the  constitutive  behavior  of  a  dry 
remolded  sand  (20-40  Ottawa  Sand).  The  sand  was  molded  into  a  thin  disk¬ 
shaped  specimen  at  the  bottom  of  a  rigid  cylindrical  test  chamber  which 
provided  lateral  constraint  (uniaxial  strain).  By  means  of  rather  ela¬ 
borate  ram  and  explosive  loading  devices,  several  specimens  were  pressure 
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loaded  with  different  rise  times  ranging  from  approximately  0.2  to  20,000 
milliseconds.  Data  for  each  test  included  the  time  history  of  the  pressure 
loading  and  the  corresponding  strain  history. 

Due  to  the  small  thickness  of  the  disk-shaped  specimen  (1.27  cm) 
intertial  effects  are  negligible  even  for  the  most  rapid  loading  rate.  That 
is,  in  reference  to  the  so-called  "multiple-reflection  theory",  the  sand- 
specimen  thickness  is  designed  to  be  sufficiently  small  to  permit  a  stress 
wave  to  multiply  propagate  back  and  forth  between  the  rigid-bottom  boundary 
and  the  free-surface  boundary  during  the  loading  rise  time.  According  to 
the  theory,  inertial  stresses  are  negligible,  and  therefore,  the  resulting 
stress-strain  histories  provide  a  direct  representation  of  the  constitutive 
behavior.  This  theory  was  independently  verified  by  the  WES  investigators 
for  their  test  specimens  by  a  simple  dynamic  analysis  (i.e.,  a  one¬ 
dimensional  wave  propagation  computer  program  using  the  actual  loading 
histories  and  piece-wise  linear  stress-strain  relations  determined  from 
static  tests). 

Figure  14  shows  the  pressure  loading  history  (on  stress)  along  with 
the  measured  strain  history  for  the  "slow"  loading  rate  which  has  a  rise 
time  of  15,000  mi  1 1 i seconds.  At  the  other  extreme,  Figure  15  shows  the 
stress  and  strain  histories  for  the  "rapid"  loading  rate  which  has  a  rise 
time  of  0.2  milliseconds.  The  resulting  stress-strain  curves  for  both 
loading  rates  are  shown  in  Figure  16.  Here  it  is  plainly  evident  that  the 
sand  specimen  exhibits  rate-dependent  stress-strain  behavior. 

Intermediate  loading  rates  with  rise  times  on  the  order  of  100  milli¬ 
seconds  gave  results  almost  identical  to  the  slow  loading  rate  experiment. 
This  leads  to  two  Important  observations,  (1)  the  non-linear  stress-strain 
relationship  for  the  slow  loading  rate  is  not  time  dependent  and  may  be 
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Figure  14.  Slow  loading  stress  and  strain  histories. 
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Figure  15.  Rapid  loading  stress  and  strain  histories 
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Figure  16.  Stress-strain  response  for  slow  and  rapid  loading. 


assumed  inviscid,  and  (?)  rate  effect  only  become  significant  when  the 
rise  time  approaches  the  sub-millisecond  range.  These  observations  will  be 
used  to  great  advantage  in  the  subsequent  parameter  identification  process. 

The  overall  strategy  for  determining  the  model's  parameters  (listed  in 
Figure  16)  is  based  on  the  observation  that  the  experimental  stress-strain 
curve  for  the  slow  loading  rate  is  an  in* i sc id-plast ic  response. 
Accordingly,  all  the  elastic  and  plastic  parameters  can  be  determined  from 
the  slow  loading-rate  test.  Once  these  parameters  nave  been  identified, 
the  viscous  parameters  (,  and  N)  can  be  determined  from  the  rapid  loading 
test.  That  is,  for  the  slow-loading  trial  simulations,  ■  is  taken  suf¬ 
ficiently  large  to  ensure  complete  viscoplastic  flow  (i.e.  mv  i  s<\  id-p  1  ast  u 
response).  In  other  words,  there  is  some  lower  limit  or  •  su«.h  tnat  any 
value  greater  than  this  limit  produces  identical  results.  In.-  ‘ina' 
for  ’  is  directed  by  the  r.  pid-1  oading  test,  i.e,  •  is  chose"  5,  tna  arc 
error  to  achieve  reasonable  agreement  between  the  predicted  a"d  meas .red 
stress-strain  slopes  from  rapid  loading. 

With  the  above  understanding,  the  following  parameter  1  dent • f 1  cat  1 or 
procedure  employs  the  strain  loading  history  from  liyure  14  as  input  irtu 
the  VPDRVR  program.  Identification  begins  by  selecting  elastic  parameters 
to  match  the  initial  unloading  slope  of  the  slow-loading  test,  lhis  slope 
is  an  elastic  confined  modulus  graphically  measured  as  57,000  MPa.  ;,sing 
this  value  along  with  an  assumed  value  of  Poisson's  ratio  -  .3,  the  bulk 
and  shear  modulus  are  set  once  and  for  all  as  recorded  in  Figure  16. 

Selection  of  the  failure  surface  parameters,  A,  B,  and  C  (Sandler 
form)  are  guided  by  the  observation  that  the  unloading  curve  begins  to 
exhibit  a  nonlinear  response  after  a  st-ess  reduction  between  5  to  ir  MPa. 
This  suggests  that  the  elastic  unloading  space  is  rather  small,  hence,  the 
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maximum  failure  Surface  height  A  is  Srt  a  »oftin#lly  VM  '<  *alu*  of  5 
MPa,  and  the  curvature  parameter  8  's  dete'^mea  by  trial  and  errur  tv  pro¬ 
duce  a  re-entry  point  on  tne  failure  surfa*.**  »r  ch  approximately  matches 
the  break  in  the  unloading  curve.  The  parameter  value  for  L  is  chosen 
slightly  ’ess  than  A  in  order  tnat  tne  difference.  A-C,  provides  a  ve»y 
Snail  value  for  tne  initial  elastic  spa*.'*  I'Mor  tv  loading.  This  is  in 
uon  f  opua  i  e  with  the  oDserva'ij"  tnat  "’i  '"itijl  elastic  response  is 
observed  upon  ini*ia*  luadtn,.  4-.  a  side  *. uvrv r t .  it  »as  observed  through 
numerical  exper imentat  ••>»*  that  a  non*  pronounced  bre*»  *r  tne  uft'oadlnj 
curve  tan  be  achieved  b,  reducing  A  and  increasing  8. 

cor  the  initial  ca;»  Surface  location,  t^.  is  set  at  an  arbitrary  small 
value  to  limit  the  si/e  of  the  initial  elastic  space  in  cor*ormance  with 
the  small  A-C  value  discussed  above.  The  cap  shape  parameter  A  is 
arbitrarily  set  at  2. t,.  Numerical  exper imentat ion  indicated  that  changes 
in  3  has  little  effect  on  the  a«>a)  stress-strain  Curve.  JtS  primary 
influence  is  to  increase  the  magnitude  of  lateral  stresses  as  R  increases 
(lateral  stresses  were  not  measured  in  the  experiments) . 

The  most  important  parameters  for  capturing  the  shape  of  the  stress- 
strain  loading  curve  are  the  cap  hardening  parameters  W  and  0.  Since  w 
represents  the  maximum  volumetric  viscoplastic  strain  that  can  be  achieved, 
it  was  initially  estimated  as  0.03  which  is  approximately  the  maximum  volu¬ 
metric  total  strain  observed  in  the  test.  As  previously  discussed.  It  is 
generally  recommended  to  choose  0 j X  (  •  0.5.  Decreasing  D  or  w  results  in 
increased  stress  magnitude,  i.e.,  a  steeper  stress-strain  loading  curve. 
After  several  trials,  the  final  values  selected  are  M  *  0.027  and  DIX  1 
0.047.  This  completes  the  identification  of  elastic  and  plastic  parame¬ 
ters. 
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During  the  above  identification  process,  tr,e  fluidity  parameter  t  was 
set  at  a  relatively  large  value  to  ensure  inviseiu  responses  for  the  slow 
loading-rate  test  (  •  -  1.0  »  id*4  millisecond*!).  After  the  final  elastic 
and  plastic  parameters  were  chosen,  t  was  repeatedly  reduced  to  ascertain 
at  what  value  of  i  the  slow  loading  -gins  to  exhibit  a  sma ’ 1  viscous 

response  which  did  not  differ  from  the  inviscid  response  by  more  thar  IX. 
'his  value  of  ■  was  determi  i  _»d  to  be  t.U?  *  io"4  milliseconds*'. 

Therefore,  tms  value  is  a  lowerbound  on  the  final  value  of  ■  that  may  be 
selected  to  be: t  fit  the  rapid- 1 oadi ng  test. 

The  last  step  in  the  identification  process  is  to  simulate  the  rapid- 
loading  test  using  the  strain  history  data  in  figure  IS  for  input  into  the 
VPDRVft  program  and  ascertaining  the  viscous  parameters  and  N  to  best 
match  the  rapid-1 oadi ng  stress-strain  curve,  all  other  parameters  remaining 
the  same.  Here  N  was  set  to  1.0  (not  varied)  and  the  final  choice  for  Y  is 
0.2  x  10"4  millisecond*!,  an  order  of  magnitude  greater  than  the  lower 
bound  established  above.  As  a  final  check,  the  slow-loading  test  was  rerun 
with  the  final  parameters  and  identical  results  were  obtained.  Moreover, 
l ntennedi ate- 1 oadi ng  rates  with  rise  times  on  the  order  of  100  mi  1 1 i seconds 
were  run,  and  the  resulting  stress-strain  responses  did  not  differ  signifi¬ 
cantly  from  the  slow-loading  rate  (i.e.,  in  conformance  with  experimental 
observat ions ) . 

'Jpon  examining  Figure  16,  it  is  evident  that  the  viscoplastic  cap 
model  accurately  reflects  the  test  data.  For  the  slow-loading  test,  the 
cap  surface  continually  moves  with  the  stress  state  producing  a  stiffening 
stress-strain  response  as  shown.  Since  the  stress  state  is  on  the  cap  sur¬ 
face  during  loading,  the  immediate  unloading  response  is  initially  elastic 
pnor  to  re-entering  the  failure  surface.  On  the  other  hand,  for  the  fast- 
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loading  rate,  the  cap  surface  lags  behind  the  loading  stress  state  pro¬ 
ducing  an  "apparent"  softemnig  stress-strain  response  (of  course,  this  is 
really  a  time-dependent  effect).  Just  prior  to  unloading, the  stress  state 
is  well  above  the  cap  surface,  thus  when  unloading  occurs,  the  stress  state 
remains  in  the  viscoplastic  domain  producing  additional  strain  accumula¬ 
tions  as  shown.  The  correlation  between  the  model's  performance  and  the 
observed  performance  is  truly  quite  remarkable,  particularly  with  regard  to 
matching  the  rapid-loading  behavior  characteristic  of  ground  shock 
problems. 

Summary  and  Recommendations 

Summary  and  Conclusions.  In  Part  I  of  this  report  a  theoretical  for¬ 
mulation  for  viscoplastic  tension  cutoff  was  developed  based  on  a  Jj  stress 
criterion.  For  completeness,  this  formulation  was  presented  with  the  pre¬ 
vious  CAP75  viscoplastic  formulation  (1)  providing  a  complete  description 
of  viscoplastic  behavior  for  tension  cut-off,  failure  surface,  and  cap  har¬ 
dening.  A  numerical  solution  strategy  for  the  complete  model  was  presented 
and  coded  in  the  computer  program  VPDRVR  (Appendix).  This  algorithm 
employs  a  variable  Crank-Nicolson  time  integration  scheme  together  with 
Newton-Raphson  iteration  procedure  to  solve  for  the  six-component  stress 
history  resulting  from  an  arbitrary  six-component  strain  loading  schedule. 
Also,  the  program  solves  the  inverse  problem,  i.e.,  stress  loading  input 
strain  history  output. 

The  new  tension-cutoff  algorithm  was  tested  against  an  exact  solution 
for  the  case  of  uniaxial-stepped-strain  loading.  Perfect  agreement  was 
obtained.  It  was  concluded  that  the  fluidity  parameters  in  the  tension 
domain  should  be  at  least  an  order  of  magnitude  larger  than  that  in  the 


65 


v i scop 1  a s t i c  cap  model  is  well  suited  for  capturing  the  time-dependent 
behavior  of  soils  and  rocks  over  a  wide  range  of  loadings.  Future  enhan¬ 
cements  of  the  model  can  easily  overcome  the  shortcomings  noted  above. 
Future  Recommendations .  Recommendations  for  future  efforts  are  divided 

into  two  main  areas;  "model  enhancement"  and  "automated  parameter 
identification".  With  regard  to  model  enhancement,  two  improvements  are 
suggested.  First  and  foremost  it  is  recommended  to  generate  the 
appropriate  functional  forms  of  the  model  to  provide  the  capabaility  of 
simulating  tertiary  creep.  This  could  be  done  by  introducing  a  history 
dependent  function  for  the  fluidity  parameter  and/or  a  strain  softening 
function  for  the  cap.  Sufficient  experimental  data  currently  exists  to 
meaningfully  undertake  this  enhancement.  The  second  enhancement  is  con¬ 
cerned  with  simulating  tension  damage  accumulation  associated  witn  cyclic 
loading.  Again,  this  could  be  done  with  special  functional  forms  for  the 
tension  fluidity  parameter  and/or  softening  functions  for  the  tension 
failure  surface.  However,  to  meaningfully  undertake  this  effort,  addi¬ 
tional  experimental  data  is  required. 

Lastly  with  regard  to  automated  parameter  identification,  it  is  recom¬ 
mended  to  re-structure  the  VPDRVR  program  into  an  interactive,  user- 
friendly,  identification  program.  For  "ideal"  data  the  program  would 
determine  all  the  model  parameters  with  very  little  assistance  from  user. 
For  "non-ideal"  data  a  close  interaction  between  the  user  and  the  computer 
is  the  best  approach.  Here  it  is  envisioned  that  the  user  would  specify 
several  constraints  (e.g.,  slope  of  unloading  curve,  initial  size  of 
elastic  domain,  etc.)  along  with  both  stress  and  strain  response  histories. 
A  first  estimate  of  the  parameters  would  be  determined  by  the  program  with 
an  over-ride  option  by  the  user.  Thereafter,  the  user  would  specify  one  or 
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more  parameters  to  be  optimized,  and  the  program  would  respond  with  the 
current  optimum  value  of  the  varied  parameter  along  with  diagnostic  data 
and  graphs  illustrating  the  effects  of  the  parameter.  Stepping  along  in 
this  interactive  fashion,  l.e.  changing  one  parameter  at  a  time,  a  final 
solution  can  be  obtained  in  a  matter  of  a  few  minutes,  instead  of  weeks  by 
a  batch  oriented  trial  and  error  approach.  Moreover,  the  intermediate 
diagnostic  data  is  of  tremendous  educational  value  with  regard  to 
understanding  the  model's  behavior. 

In  closing,  it  is  worthwhile  to  repeat  that  the  viscoplastic  cap  model 
has  been  shown  to  perform  extraordinarily  well  with  experimental  data  over 
a  wide  range  of  loading  environments,  as  well  as,  for  a  variety  of  geologi¬ 
cal  materials.  No  other  time-dependent  constitutive  model  has  exhibited 
this  degree  of  generality.  Accordingly,  it  is  highly  recommended  to  pursue 
the  future  development  of  this  model  along  the  lines  suggested  above. 
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APPENDIX  A 


PROGRAM  VPDRVR :  INPUT  INSTRUCTIONS 

This  Appendix  provides  input  instructions  for  the  VPDRVR  program  which 
exercises  the  viscoplastic  cap  model  with  tension  cutoff.  Only  a  very 
minor  change  to  the  original  input  instructions  (1)  are  needed  to  define 
the  tension  cutoff  parameters.  These  changes  are  on  one  card  (Group  0, 

Card  10)  which  is  extended  to  define  the  tension  fluidity  parameters  Yy  and 
Yg  and  the  hydrostatic  tension  cutoff  value  T. 

For  convenience,  the  entire  set  of  input  instructions  along  with  ten¬ 
sion  cutoff  input  is  given  here.  The  program  documentation  and  benchmark 
problems  given  in  Reference  (1)  remain  valid  and  are  not  repeated  here. 
Benchmarks  for  tension  cutoff  are  given  in  Part  I  of  this  report. 

Input  data  cards  are  grouped  in  the  following  categories: 

A.  (Cards  1  and  2):  Heading  and  Master  Control 

B.  (Cards  3,  4,  and  5):  Elastic  functions/parameters 

C.  (Cards  6,  7,  8,  and  9):  Plastic  function/parameters 

D.  (Card  10):  Viscous  functions/parameters  and  tension  cutoff 

E.  (Cards  11,  12):  Loading  schedules  for  stress  or  strain. 
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USER  INPUT  INSTRUCTIONS 


A.  Problem  Initiation,  Heading  and  Master  Control  Cards. 
Card  1.  ( 15A4)  Heading 


Columns 

Variable 

Entry  Description 

Notes 

01-60 
( 15A4) 

TITLE 

Descriptive  problem  title,  (program 
terminates  if  TITLE(l)  =  STOP). 

(1) 

Card  2. 

(415,  Al, 

2F10.0)  Master  Controls 

Columns 

Variable 

Entry  Description 

Notes 

01-05 

(15) 

LTYPE 

Loading  type  identification; 

=  0,  strain  loading. 

=  1,  stress  loading. 

(2) 

06-10 

(15) 

NTSEG 

Number  of  time  segments  to  define 
loading,  (Default  =  1,  Maximum  =  30). 

(3) 

11-15 

(15) 

ITMAX 

Number  of  Newton-Raphson  iterations, 

(Default  =  10). 

(4) 

16-20 

(15) 

KPR1NT 

Output  print  control; 

=  0,  standard  response  output 
=  1,  above  plus  iteration  parameters 
=  2,  above  plus  yield  function  values. 

=  3,  above  plus  iterative  correction  vector 
=  4,  above  plus  Jacobian  matrix. 

(5) 

20-21 

(Al) 

IPLOT 

Plot  control  for  response  data  written 
to  unit  11: 

=  Y  (YES)  Data  written  to  unit  11 
=  N  (NO)  Not  written 

(6) 

22-31 

(F10.0) 

THETA 

Crank  Nicolson  integration  parameter; 

°1  °1 

(7) 

32-41 

(F10.0) 

CONVRG 

Convergence  tolerance  for  Newton-Raphson 
iteration,  (Default  =  0.01,  i.e.  1% 
relative  error). 

(8) 
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B.  Elastic  Function  and  Parameter  Cards 


Card  3. 

Columns 

01-05 

(15) 


06-10 

(15) 


Card  4. 

Columns 

01-10 

(F10.0) 

11-20 

(F10.0) 

21-30 

(F 10.0) 


Card  5. 

Col umns 

01-10 

(F10.0) 

11-20 

(F 10.0) 

21-30 

JF10.0) 


(215)  Selection  of  Elastic  Functions 

Variable  Entry  Description  Notes 

IFBMOD  Selection  of  bulk  modulus  function,  K(Ji):  (9) 

=  1,  K( Jj )  »  BOATA(I)  ,  (linear) 

=  2,  K( J2)  =  BDATA(1)/(1-BDATA(2))‘ 

( 1-BDATA(2)*EXP(BDATA(3)*J1) ) 

(Default  =  1) 

IFSMOD  Selection  of  shear  modulus  function,  G(J2) :  (10) 

=  1,  G(J2)  =  SDATA(l),  (linear) 

=  2,  G(J2)  =  SDATA( 1)/(1-SDATA(2) )* 

( 1-SDATA(2)*EXP(-SDATA(3)*J2) ) . 

(Default  =  1) 


(7F10.0)  Bulk  modulus  parameters,  BDATA 


Variable 

Entry  Description 

Notes 

BDATA(l) 

First  bulk  modulus  parameter. 

(ID 

BDATA( 2) 

Second  bulk  modulus  parameter. 

BDATA( 3) 

Third  bulk  modulus  parameter. 

(7F10.0) 

Shear  modulus  parameters,  SDATA 

Variable 

Entry  Description 

Notes 

SDATA(l) 

First  shear  modulus  parameter 

(12) 

SDATA(2) 

Second  shear  modulus  parameter 

SDATA(3) 

Third  shear  modulus  parameter 

C.  Plastic  Function  and  Parameter  Cards 


Card  6. 

(415,  G10.0) 

Selection  of  CAP75  functions 

Columns 

Variable 

Entry  Description 

Notes 

01-05 

(15) 

1FFAIL 

Selection  of  failure  surface  function 
fp  =  /j?  +  9Fj  (Ji) : 

(13) 

=  1,  =  -FDTATA(l)  +  FDATA(2)*J1 . 

=  2,  gF.  =  -FDATA(l)  +  FDATA(2)* 

1  EXP(FDATA( 3)*J1 ) . 

(Default  =  1) 

06-10 

(15) 

IFCAPR 

Selection  of  cap  surface  ellipse  ratio  R: 

=  0,  No  cap,  just  failure  surface. 

■  1,  R  =  CDATA(l). 

=  2,  R  =  CDATA( 1 ) / ( 1  +  CDATA(2) )* 

(1.0  +  CDATA( 2)*EXP(CDATA( 3)*EL) ) . 

(14) 

11-15 

(15) 

IFHARD 

Control  of  cap  hardening: 

=  0,  No  hardening,  stationary  cap. 

=  1,  CAP75  hardening  function  is  used: 
e  =  W*(EXP(D*X)  -  1). 

W  =  HDATA(l) 

D  =  HDATA(2) 

(15) 

16-20 

(15) 

KAPTYP 

Selection  for  soil  or  rock  hardening  laws: 

=  0,  soil  material . 

=  1,  rock  material . 

(16) 

21-30 

(G10.0) 

XINITL 

Initial  location  of  cap  X  on  axis. 

(17) 

Card  7. 

(F10.0)  Failure  Surface  Parameters,  FDATA. 

Columns 

Variable 

Entry  Description 

Notes 

01-10 

(F10.0) 

FDATA(l) 

First  failure  surface  parameter. 

(18) 

11-20 

(F10.0) 

FDATA( 2) 

Second  failure  surface  parameter. 

21-30 

(F10.0) 

FDATA(3) 

Third  failure  surface  parameter. 
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*  Card 

8.  (7F10.0) 

Cap  Surface  Parameters  for  R,  CDATA 

Columns 

Variable 

Entry  Description 

Notes 

01-10 

(F10.0) 

CDATA(l) 

First  cap  R  parameter. 

(19) 

11-20 

(F10.0) 

CDATA(2) 

Second  cap  R  parameter. 

21-30 

(F10.0) 

CDATA(3) 

Third  cap  R  parameter. 

*  Card 

9.  ( 7F10.0) 

Hardening  cap  parameters,  HDATA. 

Columns 

Vari able 

Entry  Description 

Notes 

01-10 

(F10.0) 

HDATA(l) 

First  hardening  parameter,  W. 

(20) 

11-20 

(F10.0) 

HDATA(2) 

Second  hardening  parameter,  D. 

*Skip  Cards  8  and  9  if 

IFCAPR  =  0. 

D.  Viscous  Function  and  Tension-Cutoff  Parameters 

Card 

10.  (15,  6F10.0)  Selection  of  viscous  function/parameters 

Columns 

Variable 

Entry  Description 

Notes 

01-05 

(15) 

IF V I  SC 

Selection  of  viscous  function  : 

=  1,4  =  (f/AN0RM)**EXPN. 

=  2,4  =  EXP((f/AN0RM)**EXPN)  -  1. 
(Default  =  1) 

(21) 

06-15 

(F10.0) 

EXPN 

Exponent  in  0  function, 

(Default  =  1.0). 

(22) 

16-25 
(F 10.0) 

GAMMA 

Fluidity  parameter, Y  . 

(23) 

26-35 

(F10.0) 

ANORM 

Normalizing  constant  in  <P  function, 
(Default  =  max(FDATA(l),  0.01) 

(24) 

-  Card  10  continued  on  next  page 
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D.  Card  10  continued  (tension  cutoff  parameters) 


Columns 

Variable 

Entry  Description 

Notes 

36-45 

(F10.0) 

GAMMAB 

Fluidity  parameter 
tension  cutoff,  Yy 
Default  =  Y 

for 

volumetric 

(a) 

46-55 

(F10.0) 

GAMMAG 

Fluidity  parameter 
tension  cutoff,  Tq 
Default  =  10.0*»T 

for 

deviatoric 

(b) 

56-65 

(F10.0) 

TCUT 

Hydrostatic  tension 

cutoff  limit,  T 

(c) 

Notes  a. 

b,  and  c 

for  Card  10. 

(a)  To  simulate  rapid  volumetric  stress  release,  GAMMAB  (Yy)  should  be 
taken  significantly  greater  than  GAMMA  (y)  which  controls  the 
viscoplastic  flow  in  the  cap/failure  regions. 

(b)  In  order  to  have  deviatoric  stresses  release  at  the  same  rate  as  volu¬ 
metric  stresses,  set  Yq  =  9K0  Yy/G0,  where  K0  and  G0  are  bulk  and  shear 
elastic  moduli.  Typically  Yg  should  be  an  order  of  magnitude  greater 
than  Yy. 

(c)  The  tension  cutoff  value  TCUT  (or  T)  triggers  tension  cutoff  whenever 
Jj  -  T  >  0.  Accordingly,  a  sensible  choice  for  T  is  in  the  range  0^ 

T  ■  FCUT.  FCUT  is  where  the  failure  surface  intersects  the  Jj  axis. 

If  it  is  desired  to  deactivate  the  tension  cutoff  procedure  entirely, 
set  T  »  FCUT.  To  insure  that  TCUT  is  not  specified  within  the  cap 
domain,  the  program  checks  that  TCUT  L0.  If  this  is  not  satisfied, 
the  program  resets  TCUT  =  L0  and  is  noted  on  the  printed  output. 


74 


E.  Input  Loading  Schedule  and  Time  Steps. 


Repeat  card  set  11  and  12  NTSEG  times;  NS  ■  1.  NTSEG 

Card  11  (F10.0,  215)  Time  segment,  number  of  steps,  print  control. 


Columns 

Variable 

Entry  Description 

Notes 

01-10 

(F10.0) 

TS(NS) 

Time  at  end  of  segment  NS. 

(25) 

11-15 

(15) 

NUMDT(NS) 

Number  of  times  steps  within  time 
segment  NS. 

(Default  =  10) 

(26) 

16-20 

(15) 

IPRNT(NS) 

Print  interval  for  standard  output: 

=  1,  every  time  step  prints  output. 

=  n,  every  nth  step  prints. 

(Default  =  1) 

(27) 

Card  12. 

(6F10.0)  Stress  or  strain  load  vector  at  time  TS(NS). 

Columns 

Variable 

Entry  Description 

Notes 

01-10 

(F10.0) 

PLOAD(l.NS) 

an  (or  qj)  at  TS(NS). 

(28) 

11-20 

(F10.0) 

PL0AD(2,NS) 

022  (or  *-22)  at  TS(NS). 

21-30 

(F10.0) 

PL0AD(3,NS) 

033  (or  C33)  dt  TS(NS). 

31-40 

(F10.0) 

PL0AD(4,NS) 

012  (°r  ^12)  at  ts(ns). 

41-50 

(F10.0) 

PL0AD( 5, NS) 

013  (or  u13)  a<-  TS(NS). 

51-60 
( F10.0) 

PL0AD(6,NS) 

o23  (or  £  23)  at  TS(NS). 

***END  OF  INPUT  FOR  ONE  PROBLEM*** 


i 
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Commentary  Notes  with  Input  Instructions: 

1.  Problems  may  be  run  back-to-back .  Terminate  the  last  problem  by 
writing  STOP  in  columns  1  to  4. 

2.  Strain  loading  implies  the  six  components  of  strain  will  be  specified 
individually  during  the  loading  schedule.  Similarly,  stress  loading 
implies  the  six  components  of  stress  will  be  individually  specified. 

3.  For  either  stress  or  strain  loading,  NTSEG  is  the  desired  number  of 
time  segments  to  define  the  loading  histories  in  a  piecewise  linear 
fashion. 

4.  Generally  10  iterations  is  more  than  sufficient  to  achieve  convergence. 
If  convergence  is  not  achieved,  it  is  a  strong  indication  that  the  time 
step  is  too  large.  Note  that  convergence  of  the  Newton-Raphson 
procedure  does  not  guarantee  accuracy.  Accuracy  can  only  be  assured  by 
repeatable  solutions  with  smaller  time  steps. 

5.  Standard  output  includes  stress  or  strain  responses,  cap  location, 
number  iterations  to  converge,  stress  invariants,  and  type  of  response. 
For  KPRINT  >  0,  additional  information  is  given  primarily  for  debugging 
purposes. 

6.  Standard  response  data  is  written  to  unit  11  for  subsequent  plotting  on 
a  CALCOMP  plotter.  Subroutine  GRAPH  is  used  for  plotting  and  may  be 
removed  or  replaced  if  desired. 

7.  For  THETA  =  0.,  the  solution  algorithm  is  explicit  resulting  in  linear 

equations  {i.e.  no  Newton-Raphson  iteration).  For  THETA  ■*  0,  the 
algorithm  is  implicit  and  generally  more  accurate  for  a  given  time  step 
size,  but  requires  Newton-Raphson  iteration,  for  THETA  0.5,  the 
algorithm  is  unconditionally  stable.  ~ 

8.  The  convergence  tolerance,  C0NVRG,  is  tested  against  the  ratio  formed 
by  the  norm  of  the  correction  vector  for  stress  (or  strain)  divided  by 
the  norm  of  the  stress  (or  strain)  vector.  Norms  are  Euclidean. 

9.  The  nonlinear  bulk  modulus  function  given  by  IFBMOD  =  2  is  taken  from 
CAPCRIVER  (NCEL  Program).  It  is  a  function  of  Jj  (first  stress 
invariant)  and  is  treated  the  same  for  loading  or  unloading. 

Additional  functions  may  be  added  to  program  in  FUNCTION  DI(I,J). 

10.  The  nonlinear  shear  modulus  function  given  by  IFSM00  *  2  is  a  function 
of  j£,  second  deviator  stress  invariant  (see  Note  9). 

11.  For  future  program  expansion,  80ATA  is  dimensioned  to  7  to  allow  incor- 
poratin  of  higher  order  nonlinear  functions. 

12.  SDATA  is  dimensioned  to  7  (see  above). 

13.  For  1FFAIL  =  1,  the  failure  surface  is  standard  Orucker-Prager  (or  Von 
Mices  if  FDATA(2)  =  0.0).  For  IFFAIL  =  2,  the  failure  surfa'-"'  is  the 
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exponential  form  suggested  by  Sandler  tor  CAWS.  Additona)  functional 
forms  may  be  added  to  the  program  in  FUNCTION  FG1. 


14.  By  setting  IFCAPR  *  0,  the  plasticity  model  is  governed  by  only  the 
failure  surface.  For  IFCAPR  *  1  or  2  the  cap  surface  is  included  with 
R  given  by  the  corresponding  functional  form.  Additional  functional 
forms  for  R  may  be  added  to  program  in  FUNCTION  FRCAP.  (Note  for 
IFCAPR  *  2,  R  -  R(EL)  where  El  is  “L*  of  cap). 

15.  If  desired,  a  nonhardening  cap  surface  may  be  used  by  setting  IFHAKD  * 

0.  Otherwise  the  CAP75  hardening  function  is  employed.  New  hardening 
functions  can  be  employed  by  modifying  SUBROUTINE  CAP76. 

16.  See  Part  1  for  the  special  hardening  rules  for  soils  (KAPTYP  «  0). 

17.  The  initial  location  of  X  defines  the  starting  position  of  the  «.ap  sur¬ 
face.  The  program  checks  that  XINITl  is  not  greater  than  FCUT,  i.e. 
the  intersection  of  the  failure  surface  with  Jl  axis.  If  it  is,  XINITl 
is  automatically  reset  slightly  less  than  FCUT.  Note,  the  so-called 
Von  Mises  Transition  employed  by  Sandler  is  not  included  in  this  devel¬ 
opment.  Thus,  if  it  is  desired  to  obtain  steady-state  viscoplastic 
solutions  to  exactly  match  CAP75  plasticity  solutions,  XINITL  should  be 
chosen  so  that  the  initial  l  location  is  not  greater  than  zero. 

18.  The  "standard  Sandler"  CAP75  failure  surface  is  the  form  given  by 
IFFAIL  =  2.  In  which  case  FDATA(l)  =  A,  FDATA( 2)  =  C,  and  FDf TA(3y  • 

B. 

19.  The  "standard  Sandler"  CAP75  cap  surface  parameter  is  the  form  gwen  by 
IFCAPR  =  1,  i.e.,  CDATA(l)  =  R. 

20.  If  IFHARD  =  0,  HDATA(l)  and  HDATA(2)  are  read  but  not  used.  If  IFCAPR  = 

0,  cards  8  and  9  are  not  read.  HDATA  as  well  as  FDATA  and  CDATA  are 

dimensioned  to  7  for  future  program  expansion. 

21.  For  geological  materials  IFVISC  =  1  is  general  i.,  the  most  popular  furm 
for  the  viscous  function.  Additional  functional  forms  may  be  added  to 
the  program  in  SUBROUTINE  PH  IF. 

22.  EXPN  need  not  be  a  whole  number,  but  must  be  greater  than  zero. 

23.  GAMMA  has  units  of  inverse  time,  the  units  (e.g.  seconds,  hours,  years) 

correspond  to  the  loading  time  units  TS  in  Card  11. 

24.  Generally  the  default  value  fo  AN0RM  is  appropriate  providing  FDATA(l)  t 

0.0  .  AN0RM  should  not  be  viewed  as  an  independent  material  parameter 
since  it  is  always  associatd  with  GAMMA  in  the  quotient  GAFf1A/AN0RM**EXPN. 

25.  Up  to  30  time  segments  may  be  used  to  define  a  piecewise  continuous 
collection  of  straight  lines  to  define  loading.  For  the  first  time 
segment,  the  program  automatically  assumes  initial  time  is  zero,  i.e. 

TS(0)  =  0.0.  Thus,  TS(1)  is  the  time  at  the  end  of  first  segment, 

TS(2)  is  the  time  at  the  end  of  the  second  segment,  etc.  Successive 
values  of  TS(NS)  must  be  greater  than  the  previous  value. 


26.  Any  number  of  time  steps  may  be  assigned  to  each  time  segment. 
Accuracy/stability  is  controlled  by  the  time  step  size  so  that  it  is 
good  practice  to  repeat  solutions  by  doubling  the  value  of  NUMDT(NS) . 
Although  the  time  step  size  may  be  specified  differently  in  each  time 
segment,  it  is  good  practice  not  to  make  changes  in  t  between  segments 
by  a  factor  of  more  than  2. 

27.  The  printout  interval  may  be  specified  differently  for  each  time 
segment . 

28.  Loading  values  at  the  end  of  each  time  segment  are  specified  indivi¬ 
dually  for  each  vector  component  of  strain  if  LTYPE  =  0,  or  each  vector 
component  of  stress  if  LTYPf  =  1.  For  the  first  time  segment  the 
initial  loading  and  responses  are  automat ical ly  assumed  zero  i.e., 

J ( 0 }  =  e(0)  =  0.  Standard  continuum  mechanics  sign  conventions  are 
observed  for  all  input  and  output.  For  example  if  a  uniaxial  stress 
loading  cycle  is  desired  in  which  is  compressed  at  a  constant  rate 
to  a  stress  value  -10.0,  held  constant,  then  reverse  loaded  at  a 
constant  rate  to  a  tensile  stress  value  of  +1.0,  and  again  held 
constant;  we  infer  NTSEG  =  4,  and  Oj]  is  described  by: 

PLOAD(l.l)  -  -10.0 
PL0AD( 1,2)  »  -10.0 
PL0AD(1,3)  =  +1.0 
PL0AD(1,4)  «  +1.0 

and  all  other  stress  components  (PL0AD)  are  zero. 
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